function [ tv, uv ] = multistep ( fun, y0, y1, y2, t0, tf, Nh )
% Multistep method for the scalar ODE in the form
% y'(t) = f(t,y(t)),  t \in (t0,tf)
% y(t0) = y0
% y(t0+h) = y1
% y(t0+2h) = y2 ; with h = (tf-th)/Nh
%
%  [ tv, uv ] = multistep ( fun, y0, y1, y2, t0, tf, Nh )
%  Inputs: fun        = function handle for f(t,y), fun = @(t,y) ...
%          y0, y1, y2 = initial values
%          t0         = initial time
%          tf         = final time
%          Nh         = number of time subintervals
%  Output: tv         = vector of time steps (1 x (Nh+1))
%          uv         = vector of approximate solution at times tv
%

tv = linspace( t0, tf, Nh+1 );
h = ( tf-t0 ) / Nh;

uv = zeros( 1, Nh+1 );
uv( 1:3 ) = [y0, y1, y2];

for n = 3 : Nh
    fn = fun( tv( n ), uv( n ) ); 
    fnm1 = fun( tv( n-1 ), uv( n-1 ) ); 
    fnm2 = fun( tv( n-2 ), uv( n-2 ) ); 
    uv( n+1 ) = uv( n ) + h * ( 23*fn - 16*fnm1 + 5*fnm2 ) / 12;
end

return