

function [X,t]=midpoint(f,x0,t0,T,N,x1)
%midpoint scheme x_{n+2}=x_n+2h*f_{n+1}
%INPUT
%f - vector fild - function handle to function y=f(x,t) returning the value
%    of vector field 
%x0 - initial value x0=x(t0)
%t0 - initial time 
%N - no of time steps computed
%x1 - approximation of the solution at t1=t0+h - the error should be of order h^2
%     i.e. |x(t0+h)-x1|=O(h*h)
%OUTPUT
%X matrix N+1 x dim of x0 - approximation of solution s t
%  the ith row of  X: X(i,:) approximates x(t0+i*h)' (transpose of the solution at time t_i=t0+i*h 
%t- vector of discrete times [t_0 .. t_i..t_N] (row) t_i=t0+i*h for h=(T-t0)/N
%EXAMPLE:
%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;");
%same but for 2D 
%x0=[1;1];sol=@(t) x0*exp(-t);x1=sol(h); 
% [X,t]=midpoint(f,x0,t0,T,N,x1); plot(t,X,t,sol(t)',";sol;");
%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;");
%pause(2);plot(t,X-sol(t)',";error;");%notice that the error is not smooth...

 M=length(x0);
 h=(T-t0)/N;
 X=zeros(N+1,M);
 t=linspace(t0,T,N+1);

 X(1,:)=x0';
 X(2,:)=x1';
 for k=3:N+1,
  fp=f(X(k-1,:)',t(k-1));
  X(k,:)=(X(k-2,:)' + 2*h*fp )';
 end%for
end%function
