


# testing AdamsB2

#order of convergence for a model initial value problem 
#dx/dt=-x x(0)-1 solution x(1)=exp(-1)

er=0;
n=40;
#vector field
function y=f1(x,t)
y=-x;
endfunction

for k=1:10,
 y=AdamsB2(@f1,0,1,exp(-1/n),n,1);
 erp=er; #previous error
 er=y(n+1)-exp(-1); #error
 printf("n=%8d h=%8.3g Error=%10.3g",n, 1/n,er);
  if (k>1) 
    printf("   ratio %g\n",erp/er);
  else printf("\n");
  endif
   n*=2;
endfor


#Testing AdamsB2 for the same ivp with constant preset h=1e-2; for large [0,T]

T=1; #T increasing
n=100; #constant h=0.01;
for k=1:3,
 [y,t]=AdamsB2(@f1,0,1,exp(-T/n),n,T);
 com=[";AdamsB2  h=",num2str(T/n)," T=",int2str(T),";"];
 plot(t,y,com,t,exp(-t),";exp(-t);");
 pause(3);
 n=n*4;T=T*4;
endfor 

#now we set the time interval to [0,15] and decrease h n=n*10

T=15; #T constant and  n increasing
n=100; 
for k=1:3,
 [y,t]=AdamsB2(@f1,0,1,exp(-T/n),n,T);
 com=[";AdamsB2  h=",num2str(T/n)," T=",int2str(T),";"];
 plot(t,y,com,t,exp(-t),";exp(-t);");
 pause(3);
 n*=10;
endfor 
n/=10;
#same h but for longer time interval - wrong again.... there are always oscilations 
 T=20; n=round((T/15)*n);
 [y,t]=AdamsB2(@f1,0,1,exp(-T/n),n,T);
 com=[";AdamsB2  h=",num2str(T/n)," T=",int2str(T),";"];
 plot(t,y,com,t,exp(-t),";exp(-t);");

