function [t,y] = imp_euler(f, lh, ic, stepsize)
%Usage: [t,y] = imp_euler(f, lh, ic, stepsize)
%Improved Euler's Method implementation for ODE y' = f(t,y). Here lh is the
%interal lh[1] <= t <= lh[2] over which the solution is computed and 
%y(0) = ic. 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.

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
sol = zeros(Nf+1,2); %Solution time and values
sol(1,1) = lh(1); %Initial time
sol(1,2) = ic; %Initial value

%March out in time using improved Euler.
for i=2:Nf
    sol(i,1) = sol(i-1,1) + stepsize;
    m0 = f(sol(i-1,1),sol(i-1,2));
    w = sol(i-1,2) + stepsize*m0;
    m1 = f(sol(i,1),w);
    m = 0.5*(m0+m1);
    sol(i,2) = sol(i-1,2) + stepsize*m;
end

%Last step, might be special.
step2 = finalT - sol(Nf,1);
sol(Nf+1,1) = finalT;
m0 = f(sol(Nf,1),sol(Nf,2));
w = sol(Nf,2) + stepsize*m0;
m1 = f(sol(Nf+1,1),w);
m = 0.5*(m0+m1);
sol(Nf+1,2) = sol(Nf,2) + step2*m;
sol(Nf+1,1) = finalT;

t = sol(:,1);
y = sol(:,2);
