import ode;
import graph;
import math;
size(200,200,IgnoreAspect);

real f(real t, real y) {return cos(y);}
//real f(real t, real y) {return 1/(1+y);}
typedef real function(real,real);

real a=0;
real b=1;
real y0=0;

real L[]={1,2};

int M=L.length; // Number of modes.

//real Y0[]=array(M,y0);
real Y0[]=new real[] {-1,2};

real[] F(real t, real[] y) {
  return sequence(new real(int m) {return f(t,y[M-m-1]);},M);
  //  return new real[] {exp((L[1]-1)*t)*y[1],
  //      -exp(-(L[1]-1)*t)*y[0]};
  //  return new real[]{-y[0]^2};
}

real[] G(real t, real[] y) {
  return F(t,y)-sequence(new real(int m) {return L[m]*y[m];},M);
}

real lambda=sqrt(0.5);
real[] tau,error,error2;
int n=25;

real order=3;

for(int i=0; i < n-1; ++i) {
  real dt=(b-a)*lambda^(n-i);
  Solution S=integrate(Y0,L,F,a,b,dt,dynamic=false,0.0002,0.0004,ERK3BS,verbose=false);
  real maxnorm=0;

  Solution E=integrate(Y0,G,a,b,1e-2*dt,dynamic=false,0.0002,0.0004,RK5);
  real[] exact=E.y[E.y.length-1];

  //  real[] exact=new real[] {exp(-b)*sin(b),exp(-L[1]*b)*cos(b)};
  for(int m=0; m < M; ++m)
    maxnorm=max(maxnorm,abs(S.y[S.y.length-1][m]-exact[m]));
  if(maxnorm != 0) {
    tau.push(dt);
    //      error.push(dt^-(order+1)*maxnorm);
    error.push(maxnorm);
  }
}

/*
  for(int i=0; i < n-1; ++i) {
  real dt=(b-a)*lambda^(n-i);
  real maxnorm=0;
  for(int m=0; m < M; ++m) {
  solution S=integrate(Y0[m],L[m],f,a,b,dt,dynamic=false,0.000,1000,RK4_375,verbose=false);
  maxnorm=max(maxnorm,abs(S.y[S.y.length-1]-exact[m]));
  }
  error2.push(dt^-order*maxnorm);
  }
*/

//scale(Log,Log);
scale(Log,Linear);

//draw(graph(tau,error),marker(scale(0.8mm)*unitcircle,red));
//draw(graph(tau,error2),marker(scale(0.8mm)*unitcircle,blue));

int[] index=sequence(error.length-1);
real[] slope=log(error[index+1]/error[index])/log(tau[index+1]/tau[index]);
real[] t=sqrt(tau[index]*tau[index+1]);
//write(t,slope);
draw(graph(t,slope),red);



xaxis("$\tau$",BottomTop,LeftTicks);
yaxis("$e/\tau^"+string(order)+"$",LeftRight,RightTicks);