
function [x,t]=AdamsBashforth19(f,t0,T,x0,x1,N=20)
%MATLAB: change the default value of  the last parameter i.e. of N (remove =20)
%implementation of 2step explicit Adams-Bashforth scheme 
%x(k+1)=x(k)+0.5*h*(3*f(k)-f(k-1)) f(k-1)=f(x(k-1),t(k-1))
%h=(T-t0)/N
%solving IVP dx/dt=f(x,t) x(t0)=x0
%INPUT
%f vector field f(x,t)
%t0 initial time
%T - end time (end of interval)
%x0 initial value i.e x(t0)=x0
%x1 approximation of x(t0+h)=x(t(1)) we should have ||x1-x(t(1))|| = O(h^2)
%N - number of compute points x(k) approx x(t0+(k-1)*h) k>1
%OUTPUT
% x=[x(1),x(k),....,x(N+1)]' approximated solution
%   kth row represents transpose of  the approximation of 
%   x(t0+(k-1)*h)
%t vector with the discrete times t=[t0 t0+h ... T]
%EXAMPLE
%a=1;f=@(x,t) a*x;N=20;t0=0;T=1.1;x0=1;x1=exp(a*(T-t0)/N);[x,t]=AdamsBashford19(f,t0,T,x0,x1,N);plot(t,x,";AB;",t,exp(a*t))
%a=-1;f=@(x,t) a*x;N=20;t0=0;T=1.1;x0=1;x1=exp(a*(T-t0)/N);[x,t]=AdamsBashford19(f,t0,T,x0,x1,N);plot(t,x,";AB;",t,exp(a*t))

  h=(T-t0)/N;
  m=length(x0);
  t=linspace(t0,T,N+1)';
  x=zeros(m,N+1);
  x(:,1)=x0;
  x(:,2)=x1;
  fk=f(x(:,1),t(1));
  for k=2:N,
     fp=fk;
     fk=f(x(:,k),t(k));
     x(:,k+1)=x(:,k)+0.5*h*(3*fk-fp);
  end   
  x=x'; %for conformity to the octave black box ode solver
end
