

1==1;
# testing order of local error (rzad schematu) of 2step Adams scheme of order 2
#it is possible ONLY if we know the exact analytic solution of IVP
#(mozliwe o ile znamy analityczne rozwiazanie)

#order of local order of the Adams scheme of order 2  for a model initial value problem (IVP)
#dx/dt=x x(0)=1 solution x(t)=exp(t)

#vector field
function y=f1(x,t)
y=x;
endfunction

sol=@(x) exp(x);

#For known solution of ODE: x'=f1(x,t) 
#we compute local error of explicit Adams cheme of order 2:
# e(t,h)=|x(t+2h)-x(t+h) -0.5h(3f(t+h,x(t+h))-f(t,x(t)]| for different t and h
#and then compute the ratio R(t;h)=e(t,h)/e(t,h/2) if
#as theory predicts e(t,h)=O(h^3) then R(t;h) should tend to 8=2^3  
tt=[0.1,1,10,100];
j=0;M=10;
err=zeros(M,length(tt));
for t=tt,
 j++;
 er=1.0; 
 h=0.5;
 for k=1:M,  
  erp=er; #previous error
  er=err(k,j)=abs(sol(t+2*h)-sol(t+h)- 0.5*h*(3*sol(t+h)-sol(t))); #local error
  printf("[%2d] t=%7.2g h=%8.3g Er=%10.3g RelE=%10.3g",k,t,h,er,er/sol(t));
  if (k>1) 
    printf("   ratio %10.4g\n",erp/er);
  else printf("\n");
  endif
   h*=0.5; #h=h/2
 endfor
endfor 
#we can compute the ratios in a compact way
format short;
disp("Ratios R(t;h)");
R=err(1:M-1,:)./err(2:M,:)
disp("many digits are the same... show the numbers with more digits");
format long e;R

#Please note that for large or small values of sol(t) we are really interested in
#the relative local error (the 4th column - RelE)

#repeat the tests but for different T or  ODE e.g. y'=-y y(0)=1 with T=1,10 etc
