function [t,x] = ode3ms(odefile,tspan,x01,nsteps)
% [t,x] = ode3ms('odefile',tspan,x0,nsteps)
% Integrates x'=f(t,x) by 2-step explicit 3rd-order formula.
%
% Usage (requires zdv.m):
% >> tspan = [0 2]; x0 = [sqrt(2) sqrt(2)]; nsteps=26;
% >> [t,x]=ode3ms('zdv',tspan,x0,nsteps);
% >> plot(t,x), grid, axis([0 2 -5 5])
% >> xlabel('t'), ylabel('x(t)')
% >> title('Solution of x''=0, x(0)=sqrt(2) by 2-step 3rd-order formula')

t0=tspan(1); tf=tspan(2);
t = linspace(t0,tf,nsteps+1);
neq = length(x01(:))/2;
x01=reshape(x01,neq,2);
x0 = x01(:,1); x1 = x01(:,2);
fm2 = feval(odefile,t(1),x0);
fm1 = feval(odefile,t(2),x1);
x = [x0 x1 zeros(neq, nsteps-1)];
for n = 3:nsteps+1
  x(:,n) = -4*x(:,n-1)+5*x(:,n-2) + 4*fm1 +2*fm2;
  fm2    =fm1;
  fm1    =feval(odefile,t(n),x(:,n));
end