function [tn,yn] = impEuler(f,t0,tend,y0,h)
%
% Solve the ODE y'=f(t,y) with improved Euler method. 
% Input: f       : Function handle to the rhs. of the ODE. 
%        t0, tend: The integration interval. 
%        y0      : Initial value   
%        h       : Stepsize
% Output: tn     :
%         yn     : The solution at the points tn.

t = t0;    % Current solution
y = y0;

tn = t0;   % Store the solution
yn = y0;

while t < tend
    if t+h > tend;
        h = tend-t;
    end
    
    % The Improved Euler step
    k1 = feval(f,t,y);
    k2 = feval(f,t+h,y+h*k1);
        
    y = y+0.5*h*(k1+k2);
    t = t+h;
    
    % Store the solution
    tn = [tn,t];
    yn = [yn,y];
end


    
    



 


    



