
#solving odes by lsode - 1st part 
#simple problems in 1 or 2d - blow up real and false
#stiff equatiosn - explicit vs implicit methods in lsode

#dx/dt=x x(0)=1
f=@(x,t) x;
t=linspace(0,1,100);
X=lsode(f,1,t);sol=exp(t)';
plot(t,X,t,sol,"r");pause(1);
er=abs(sol-X);
plot(t,er,";error;");pause(1);


#2d example linearized pendulum
f=@(x,t) [0,1;-1,0]*x;
t=linspace(0,7,1000);
X=lsode(f,[0;1],t);
plot(t,X); pause(1);

plot(X(:,1),X(:,2));pause(1);

#zobaczmy czy okrag sie zachowa
T=6500;
X=lsode(f,[0;1],[0,T]);
norm(X(:,2),2)

#blow up x'=x^2 x(0)=x0 >0 
#sol x(t)=1./(1/x0-t)
f=@(x,t) x*x;
ep=6e-7;
t=linspace(0,1-ep);
X=lsode(f,1,t);
plot(t,X); pause(1);

#false blow-up numerical
#x'=800*x x(0)=1
f=@(x,t) 700*x;
t=linspace(0,1);
X=lsode(f,1,t);
exp(700)
#if we take x'=800x we would get a  false blow up - fl overflow...
#f=@(x,t) 700*x;
#t=linspace(0,1);
#X=lsode(f,1,t)
#exp(700)

#stiff system 
A=[998 1998; -999 -1999]
eig(A) #so very stiff
f=@(x,t) A*x;
tic;lsode(f,[1;1],[0,20]);ts(1)=toc
tic;lsode(f,[1;1],[0,60]);ts(2)=toc
lsode_options("integration method","adams")
tic;lsode(f,[1;1],[0,20]);tns(1)=toc
tic;lsode(f,[1;1],[0,60]);tns(2)=toc
tic;lsode(f,[1;1],[0,200]);toc
#on [0,200] explicit Adamsschems fail
tns./ts
#but bdf not
lsode_options("integration method","stiff")
tic;lsode(f,[1;1],[0,200]);toc


