function [t,y] = impeuler_sys(f, lh, ic, stepsize)
%Usage: [t,y] = impeulers_sys(f, lh, ic, stepsize)
%Improved Euler's Method implementation for ODE system y' = f(t,y) 
%hence y can be a vector. Here lh is the interal lh[1] <= t <= lh[2] over 
%which the solution is computed and y(lh[1]) = ic, where ic is the vector
%of initial conditions at time t = lh[1]. The step size is in argument 
%"stepsize". If the number of steps from initial to final time is not an 
%integer, the last step is incremental of proper length. Solution times 
%are in t, solution values in y.

%Prep for marching out in time
finalT = lh(2);
N = floor((finalT-lh(1))/stepsize); %Round down
T1 = lh(1)+N*stepsize; %Corresponding time
Nf = N+1; %Number of steps we'll take, unless T1 = finalT
if abs(T1-finalT) < 1.0e-8 %To some tolerance
    Nf = N;
end

%Make that ic is a row vector.
ic = ic(:)';

%Set up arrays for solution return
N = length(ic);
t = zeros(Nf+1,1);
y = zeros(Nf+1,N);
t(1) = lh(1); %Initial time
y(1,:) = ic; %Initial value

%March out in time
for i=2:Nf
    t(i) = t(i-1) + stepsize; %Update time
    
    %Extrapolate solution. First one Euler step
    df = f(t(i-1),y(i-1,:));
    df = df(:)'; %Make sure df is a row vector
    w = y(i-1,:) + stepsize*df;
    %And another Euler step
    df = f(t(i),w);
    df = df(:)';
    w2 = w + stepsize*df;
    y(i,:) = 0.5*(y(i-1,:)+w2); %Average base and double extrapolated
end

%Last step, might be special/incremental.
t(Nf+1) = finalT; %Update time
df = f(t(Nf),y(Nf,:));
df = df(:)';
w = y(Nf,:) + (finalT-t(Nf))*df;
df = f(t(Nf+1),w);
df = df(:)';
w2 = w + (finalT-t(Nf))*df;
y(Nf+1,:) = 0.5*(y(Nf,:)+w2);