
1==1;
%check the order of LTE of forward Heun scheme
%x(n+1)=x(n) + 0.5*h*(f(x(n),t(n))+ f(x(n)+h*f(x(n),t(n)),t(n)+h))

function [e,r]=LTEforHeun(f,x,t,h0=0.1,M=8)
%checking the order of LTE of Heun at given t
%INPUT
%f - vector field i.e the rhs function f(x,t)
%x - solution of x'=f(x,t) 
%t time at which we compute e(t) LTE of the Heun scheme
%h0 - initial step
%M - no of halving the step
%OUTPUT
%e - the vector of LTE of all h's  for given t 
%r the vector with ration e(h)/e(h/2)
 h=h0;
 for k=1:M,
  ft= f(x(t),t);
  e(k)=norm(x(t+h)-x(t)-0.5*h*(ft+f(x(t)+h*ft,t+h)),'inf');
  h*=0.5;
 endfor
 e=e';
 r=e(1:M-1)./e(2:end);
end
lso
x=@(t) exp(a*t);
disp("LTE of Heun for y'=y");
[e,r]=LTEforHeun(f,x,t,h0=0.1,M=8)


a=-1;
f=@(x,t) a*x;
x=@(t) exp(a*t);
disp("LTE of Heun for y'=y");
[e,r]=LTEforHeun(f,x,t,h0=0.1,M=8)

% penduulum equation x''=-x 
%x(t)=[cos(t);-sin(t)]
%x1'=x2 x2'=-x1 or x'=[0,1;-1,0]*x
f=@(x,t) [0,1;-1,0]*x;
x=@(t) [cos(t);-sin(t)];
disp("LTE of Heun for y''=-y");
[e,r]=LTEforHeun(f,x,t,h0=0.1,M=8)

t=30;
disp("LTE of Heun for y''=-y");
[e,r]=LTEforHeun(f,x,t,h0=0.1,M=8)