#!/usr/bin/octave

1==0;

#the old version of AdamsB2() function !!!!!
function [X,t]=AdamsB2(pw,t0,x0,x1,h=0.1,N=100)
#explicit Adams scheme of 2nd order solves approximately
#x'=pw(x,t) with x(t0)=x0
#y=pw(x,t) fucntion handle to the vector field  - output: vector y of size of x
#t0- initial time  
#x0- vector of length M  (vertical) - initial value
#x1 approximation  of solution at  t_1=t0+h (2 step scheme)
#Output (as in lsode())
#X- matrix Nx M which has  an approximated solution at time moment t(k):
#x(t(k)) in the kth row for t(k)=t0+(k+1)*h
#t -vector with time moments t(k)
M=length(x0);
X=zeros(N,M);
X([1 2],:)=[x0'; x1'];
fx=pw(X(1,:),t0); #f(x0.t0)
t=t0+(0:N-1).* h;
for k=3:N
 pfx=fx;
 fx=pw(X(k-1,:),t); 
 X(k,:)=X(k-1,:)+0.5*h*(3*fx'-pfx');
endfor

endfunction


#testing Adams-Bash scheme of 2nd order with different staring values x1
#y'=y  and y(0)=1

f=@(x,t) x;
x0=1;

T=10;t0=0;
erre=errh=errx=0; N=20;
sol=exp(10);
printf("  h     &heun  start&e(h)/e(h/2)&exEuler start&ratio&exact x1\n\n");
for k=1:9,
 h=(T-t0)/(N-1); 
 [Xx]=AdamsB2(f,t0,x0,exp(h),h,N);
 [Xe]=AdamsB2(f,t0,x0,x0+h*f(x0,t0),h,N); #start x1 - one step explicit Euler
 x1h=x0+0.5*h*(f(x0,t0)+f(x0+h*f(x0,t0),t0));#one step of Heun 
 [Xh]=AdamsB2(f,t0,x0,x1h,h,N);
 errep=erre;errhp=errh;errxp=errx;
 errh=abs(Xh(N)-sol);#Heun
 erre=abs(Xe(N)-sol);# explicit Euler
 errx=abs(Xx(N)-sol); #exact

 printf("%2.1e &  %2.2e &  %2.2f ",h,errh,errhp/errh);
 printf("& %2.2e &  %2.2f ",erre,errep/erre);
 printf("& %2.2e &  %2.2f \\\\ \n",errx,errxp/errx);
 N*=2;
endfor



#testing Adams-Bash scheme of 2nd order with different staring values
#y'=-y and y(0)=1
disp("Same for y'=-y y(0)=1 -> y(1)");



f=@(x,t) -x;
x0=1;

T=1;t0=0;
erre=errh=errx=0; N=20; 
sol=exp(-1);
printf("  h     &Heun  start&e(h)/e(h/2)&exEuler start&ratio&exact x1\n\n");
for k=1:9,
 h=(T-t0)/(N-1); 
 [Xx]=AdamsB2(f,t0,x0,exp(-h),h,N);
 [Xe]=AdamsB2(f,t0,x0,x0+h*f(x0,t0),h,N); #start x1 - one step explicit Euler
 x1h=x0+0.5*h*(f(x0,t0)+f(x0+h*f(x0,t0),t0));#one step of Heun 
 [Xh]=AdamsB2(f,t0,x0,x1h,h,N);
 errep=erre;errhp=errh;errxp=errx;
 errh=abs(Xh(N)-sol);#Heun
 erre=abs(Xe(N)-sol);# explicit Euler
 errx=abs(Xx(N)-sol); #exact

 printf("%2.1e &  %2.2e &  %2.2f ",h,errh,errhp/errh);
 printf("& %2.2e &  %2.2f ",erre,errep/erre);
 printf("& %2.2e &  %2.2f \\\\ \n",errx,errxp/errx);
 N*=2;
endfor
