
1==1;
graphics_toolkit("gnuplot");
f=@(x,t) x; %function handle to a vector field: (x,t)|-->x
t0=0;
x0=1;
T=1;
N=10;
[x,t]=exEuler19(f,t0,T,x0,N)
plot(t,x,t,exp(t));pause(2)

%we test it for dx/dt=x but in 2D
x0=[1;1];
[x,t]=exEuler19(f,t0,T,x0,N);
f=@(x,t) diag([-1,1])*x;
x
[x,t]=exEuler19(f,t0,T,x0,N);
x


function er=testeelin19(a=1,M=6)
%tests of the order of convergence of explicit Euler for 
%y'=ay x(0)=1 on [0,T] a=1,-1,10,-10
%T=1,4,16,64

T=1;
f=@(x,t) a*x; %function handle to a vector field: (x,t)|-->a*x
t0=0;
x0=1;
for k=1:3,
 N=10;
 for l=1:M,
   [x,t]=exEuler19(f,t0,T,x0,N);
   com=[";a=",num2str(a)," T=",num2str(T)," N=",num2str(N),";"];
   plot(t,x,com,t,exp(a*t));pause(1);
   err=er(l,k)=norm(exp(a*T)-x(N+1));
   printf("a=%G T=%g N=%g er=%g\n",a,T,N,err);
   N*=2;
 end
   %er(1:M-1,k)./er(2:M,k)
   T*=4;
end 
end

%a=1
er=testeelin19(a=1,M=6)
er(1:M-1,:)./er(2:M,:)

%a=-1
er=testeelin19(a=-1,M=6)
er(1:M-1,:)./er(2:M,:)

%2D example linear penduulum
f=@(x,t) [x(2);-x(1)]
x0=[1;0];
t0=0;
T=7;
N=140;
[x,t]=exEuler19(f,t0,T,x0,N);
erT=x(end,:)-[cos(T),-sin(T)]
plot(t,x)%[cos(t),-sin(t)]
pause(1);
%trajectory (orbit?)
plot(x(:,1),x(:,2),cos(t),-sin(t));
pause(1);
%smaller h div 2
N*=2;
[x,t]=exEuler19(f,t0,T,x0,N);
erT=x(end,:)-[cos(T),-sin(T)]
plot(t,x)%[cos(t),-sin(t)]
pause(1);
%trajectory (orbit?)
plot(x(:,1),x(:,2),cos(t),-sin(t))
