
function [x,t]=adamsb2s17(f,x0,x1,t0,T,n=100)
%prosta implementacja otwartego schematu Adamsa-bashforda 2 krokowego
%explicit (forward) Adams Bashford 2 step scheme
%x(n+1)=x(n)+0.5*h*(3f(n)-f(n-1))
%INPUT
%f pole wektorowe /vector field f(x,t)
%x0 warunek poczatkowy/initial value
%x1 - starting value aproximating x(t0+h) computed by a one-step scheme in practise
%t0 start czasu/inital time
%T - koniec odcinka calkowania
%n - ilosc pktow czasowych tj krok czasowy h=T/n
%OUTPUT
% x - wektor dlugosci n+1 z przyblizeniem rozwiazania dx/dt=pw(t,x) x(t0)=x0
%     obliczonymi schematem Adamsa-Bashforda 2 krokowym
%     x(k+2)=x(k+1)+0.5*h*(3*f(t(k+1),x(k+1))-f(t(k),x(k)));
%     x(k) przybliza x(t(k)) x - rozwiazanie dokladne
%     vector of length n+1 approximating solution of IVP: dx/dt=pw(t,x) x(t0)=x0
%     COMPUTED BY 2step Adams-Basford scheme:  i.e. 
%     x(k+2)=x(k+1)+0.5*h*(3*f(t(k+1),x(k+1))-f(t(k),x(k)));
%     x(k) approximates x(t(k)) x - the solution of IVP
%t - discrete times t=[t0,t0+h,..T]
%Przyklad/Example
%n=100;t0=0;T=1;h=(T-t0)/n;pw=@(x,t) -x;sol=@(t) exp(-t);
%[x,t]=adamsb2s17(pw,x0=1,x1=sol(h),t0,T,n);plot(t,x,";AdamsB2s;",t,s=sol(-t),";exp(-t);");pause(2);plot(x-s,";error;");
 x=zeros(length(x0),n+1);
 t=linspace(t0,T,n+1)';
 h=(T-t0)/n;
 x(:,1)=x0;
 fkp1=f(x(:,1),t(1));
 x(:,2)=x1; %if x1 not given: x(:,2)=x0+h*fkp1;
 for k=1:n-1,
  fk=fkp1;
  fkp1=f(x(:,k+1),t(k+1));
  x(:,k+2)=x(:,k+1)+0.5*h*(3*fkp1-fk);
 end
 x=x'; %to stay consistent with a standard octave ODE solver lsode()
end
