
#testing if midpoint works at all
f=@(x,t) -x; 
x0=1;t0=0;T=1;N=10;
h=(T-t0)/N;
sol=@(t)exp(-t);
x1=sol(h); 
[X,t]=midpoint(f,x0,t0,T,N,x1); 
plot(t,X,t,sol(t),";sol=exp(-t);");pause(2);
%same but for 2D 
x0=[1;2];sol=@(t) x0*exp(-t);
x1=sol(h); 
 [X,t]=midpoint(f,x0,t0,T,N,x1); 
plot(t,X,t,sol(t)',";sol exp(-t);");pause(2);
%real 2d example
f=@(x,t) [0,1;-1,0]*x; x0=[1;0];t0=0;T=1;N=10;h=(T-t0)/N;sol=@(t)[cos(t);-sin(t)];
x1=sol(h);  [X,t]=midpoint(f,x0,t0,T,N,x1); 
plot(t,X,t,sol(t)',";sol=[2*exp(-t),exp(-t)];");
plot(t,X-sol(t)',";error;");pause(2);
%notice that the error is not smooth...

%testing order of convergence of midpoint for dx/dt=-x x(0)=1
f=@(x,t) -x; 
x0=1;t0=0;T=1;N=10;
h=(T-t0)/N;
sol=@(t)exp(-t);
x1=sol(h);er=0;
for k=1:10,
 h=(T-t0)/N;
 x1=sol(h);
 y=midpoint(f,x0,t0,T,N,x1); %midpoint(@f1,0,1,exp(-1/n),n,1);
 erp=er; %previous error
 er=y(N+1)-sol(T); %error
 printf("n=%8d h=%8.3g Error=%10.3g",N, h ,er);
  if (k>1) 
    printf("   ratio %g\n",erp/er);
  else printf("\n");
  end%if
   N*=2;
end%for

%testing midpoint for dx/dt=-x x(0)=1 for [0,T] constant h T longer and longer


f=@(x,t) -x; x0=1;t0=0;
N0=N=10;
h=1/N;
T=t0+N*h;
sol=@(t)exp(-t);

for k=1:7,
x1=sol(h); 
[X,t]=midpoint(f,x0,t0,T,N,x1);
com=[";T=",num2str(T)," h=",num2str(h),";"] ;
plot(t,X,com);pause(2);
N=N+N0*2; 
T=t0+N*h; %T=T+2;
end
#same for h=1/50

N0=N=50;h=1/N;
T=t0+N*h;
sol=@(t)exp(-t);

for k=1:7,
x1=sol(h); 
[X,t]=midpoint(f,x0,t0,T,N,x1);
com=[";T=",num2str(T)," h=",num2str(h),";"]; 
plot(t,X,com);pause(2);
N=N+N0*2; 
T=t0+N*h; %T=T+2;
end
%try smaller h.... always for some T we get oscillations due to 
%the lack of strong stability of midpoint... 


 
%2D example:
vf=@(x,t) [0,1;-1,0]*x;
a=0;b=8.3;
N=round((b-a)*5);x0=[1;0];
sol=@(t) [cos(t);-sin(t)];
x1=sol(a+(b-a)/N);
[x,t]=midpoint(vf,x0,a,b,N,x1);%midpoint(vf,a,x0,x1,N,b);
solt=sol(t)';
plot(t,x,";exE;",t,solt,";sol;");
pause(3);
com=[";midpoint err h=",num2str(b/N,"%.2e")," b=",int2str(b),";"];
plot(t,x-solt,com);
pause(3);
%now the convergence order by halving the step 
N=70;a=0;b=2*pi;x0=[1;0];
r(1)=1;
for k=1:7,
 NN(k)=N;
 x0=sol(a);x1=sol(a+(b-a)/N);
 [x,t]=midpoint(vf,x0,a,b,N,x1);%midpoint(vf,a,x0,x1,N,b);
 er(k)=norm(x(N+1,:)-sol(b)');
 if(k>1) r(k)=er(k-1)/er(k);end%if
 N*=2;
end%for
printf("2d order example: linear penduluum: testing order by the halving step method \n");
[(b-a)./NN;er;r;log2(abs(r))]'
%1st col: h; 2nd col: err; 3rd col: ratio (should be close to 4) and 4th col log_2 ratio i.e. order (should be close to 2)
solt=sol(t)';
plot(t,x,";midp;",t,solt,";sol;"); 
pause(3);
com=[";midpoint err h=",num2str(b/N,"%.1e")," b=",int2str(b),";"];
plot(t,x-solt,com);pause(3);
com=[";midpoint h=",num2str(b/N,"%.1e"),";"];
plot(x(:,1),x(:,2),com,solt(:,1),solt(:,2),";sol;");pause(3);
%Okay 


