

function [y,t]=AdamsB2step(vf,t0,x0,x1,n,t)
%Function solves approx. dx/dt=vf(x,t);x(t0)=x0
%using Adams-Bashford 2step scheme i.e. computes x(t0+k*h) k=0,..,n
%Output
%y=[y1 ... y(n+1)] -> y(k+1) approximate sol(t0+k*h)
%t=[t0, t0+h...t0+n*h]
%Input
%vf  - function vf(x,t)- vector field - with input x and t and ouput value of vf(x,t)
%t0, x0 values from initial value
%x1 approximation of x(t0+h) for h=(t-t0)/n
%n- number of steps 
%t - time value
%EXAMPLE:
%vf=@(x,t) x; sol=@(t) exp(t);N=10;a=0;b=1;x0=1;[x,t]=AdamsB2step(vf,a,x0,sol(a+(b-a)/N),N,b);plot(t,x,t,sol(t),";sol;")
%2D example:
%vf=@(x,t) [0,1;-1,0]*x;sol=@(t)[cos(t);-sin(t)];N=40;a=0;b=8.3;x0=[1;0];[x,t]=AdamsB2step(vf,a,x0,sol(a+(b-a)/N),N,b);
%norm(x-sol(t)','inf');plot(t,x,";AdamsB2s;",t,sol(t)',";sol;");pause(3);plot(x(:,1),x(:,2));

h=(t-t0)/n;
t=linspace(t0,t,n+1);
y=zeros(length(x0),n+1);
y(:,1)=xpp=x0;y(:,2)=xp=x1;
fp=vf(xpp,t0);
 for k=2:n,
  fpp=fp;
  fp=vf(xp,t(k)); 
  xn=xp+0.5*h*(3*fp-fpp);
  y(:,k+1)=xn;
  xpp=xp;xp=xn;
 end%for
 y=y'; %conforming to lsode
end%function

