

%solving simple IVP y'=ay y(0)=1 a=-1,1,10,-10,100,100 on [0,T]= T=1,10,100
graphics_toolkit("gnuplot");
T=1;
a=1;
t=linspace(0,T,200);
X=lsode(@(x,t) a*x,1,t);
plot(t,X);pause(1);
a=-1;
t=linspace(0,T,200);
X=lsode(@(x,t) a*x,1,t);
plot(t,X);pause(1);

%artificial blow-uplus
a=100;
T=10;
t=linspace(0,T,200);
X=lsode(@(x,t) a*x,1,t);

%linear penduulum
%y''=-y y(0)=0 y(1)=1 y(t)=sin(t)
f=@(x,t) [0,1;-1,0]*x;
t=linspace(0,2*pi,200);
X=lsode(f,x0=[0;1],t);
plot(t,X(:,1),";x(t);",t,X(:,2),";dx(t);");pause(1);
plot(X(:,1),X(:,2),";x vs x';");pause(1);
%check if the orbit periodic
norm(X(end,:)-x0')
%repeat for 100*PI
X=lsode(f,x0=[0;1],[0,100*pi]);
%check if the orbit periodic
norm(X(end,:)-x0')
%stiffness
%solving stiff system with lsode witth options for stiff and non-stiff ODEs


A=[998 1998;-999 -1999]
eig(A)
t=zeros(3,2);
f=@(x,t) A*x;
T=1;
lsode_options("integration method","stiff");
for k=1:3,
  tic;lsode(f,[1;1],[0,T]);t(k,1)=toc;
  T*=5;
end

 T=1;
lsode_options("integration method","non-stiff");
for k=1:3,
  tic;lsode(f,[1;1],[0,T]);t(k,2)=toc;
  T*=5;
end 

t
t(:,2)./t(:,1)


%y'=f(y,t) y(0)=s F(s)=y(1;s) the solution at t=1 with initial value equaling s
%plot y(s) for s\in[0,1] for f(y,t)=y*(1-y)*sin(t+y*y)
f=@(y,t) y.*(1-y).*sin(t+y.*y)
s=linspace(0,1);
F=zeros(size(s));
tic;
for k=1:100,
  X=lsode(f,s(k),[0,1]);
  F(k)=X(2);
end
toc  
plot(s,F);pause(1);
tic;X=lsode(f,s',[0,1]);;toc
F1=X(2,:);
norm(F-F1,"max")
plot(s,F1)