
1==1;
%check the order of LTE of forward Euler

function [e,r]=LTEforFEuler(f,x,t,h0=0.1,M=8)
%checking the order of LTE of forward Euler 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 forward Euler 
%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,
  e(k)=norm(x(t+h)-x(t)-h*f(x(t),t),'inf');
  h*=0.5;
 endfor
 e=e';
 r=e(1:M-1)./e(2:end);
end


%backward Euler

function [e,r]=LTEforBEuler(f,x,t,h0=0.1,M=8)
%checking the order of LTE of backward Euler 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 forward Euler 
%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,
  e(k)=norm(x(t+h)-x(t)-h*f(x(t+h),t+h),'inf');
  h*=0.5;
 endfor
 e=e';
 r=e(1:M-1)./e(2:end);
end

%y'=ay y(t)=exp(t)

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

a=1;
f=@(x,t) a*x;
x=@(t) exp(a*t);
disp("LTE of forward Euler for y'=y");
[e,r]=LTEforFEuler(f,x,t,h0=0.1,M=8)
disp("LTE of backward Euler for y'=y");
[e,r]=LTEforBEuler(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 forward Euler for y''=-y");
[e,r]=LTEforFEuler(f,x,t,h0=0.1,M=8)
disp("LTE of backward Euler for y''=y-");
[e,r]=LTEforBEuler(f,x,t,h0=0.1,M=8)

t=30;
disp("LTE of forward Euler for y''=-y");
[e,r]=LTEforFEuler(f,x,t,h0=0.1,M=8)
disp("LTE of backward Euler for y''=y-");
[e,r]=LTEforBEuler(f,x,t,h0=0.1,M=8)