

%simple ODEs (IVP)


%the simplest examples

%dx/dt=a*x x(0)=1
a=1;
T=1;
x0=1;
X=lsode(@(x,t) a*x,x0,t=linspace(0,T));
plot(t,X,t,exp(t),";exp(t);");pause(1);

%artificial blowup
a=20;
T=40;
X=lsode(@(x,t) a*x,x0,t=linspace(0,T));
plot(t,X,t,exp(a*t),";exp(at);")


%linear penduulum
%x''=-x x(0)=1 x'(0)=0 x(t)=cos(x)
T=7;x0=[1;0];
X=lsode(@(x,t) [0,1;-1,0]*x,x0,t=linspace(0,T));
plot(t,X(:,1),";x(t);",t,X(:,2),";x'(t);");
pause(1);
plot(X(:,1),X(:,2),";x vs x'(t);");pause(1);

%checking if periodicity is kept
j=1;x0=[1;0];clear er;
for k=0:500:10000,
  T=2*pi*(k+1);
  X=lsode(@(x,t) [0,1;-1,0]*x,x0,[0,T]);
  er(j,:)=[k+1,norm(X(2,:)'-x0,'inf')];
  j++;
end  
X=lsode(@(x,t) [0,1;-1,0]*x,x0,t=linspace(0,T,50*round(T)));

plot(X(:,1),X(:,2),";x vs x'(t);");pause(1);
er

%nonlinear penduulum
T=10;
X=lsode(@(x,t) [x(2);-sin(x(1))],x0,t=linspace(0,T,50*round(T)));
plot(X(:,1),X(:,2),";x vs x'(t);");pause(1);

T=30000;
X=lsode(@(x,t) [x(2);-sin(x(1))],x0,t=linspace(0,T,50*round(T)));
plot(X(:,1),X(:,2),";x vs x'(t) nonlin pend ling time;");pause(1);


%backward integration of IVP for t<t0

function X=backivp(F,x0,t)
%F vector field Y=F(x,t)
%x0 - initial value at t=t0
%t- backward times vector t=[t(1),..,t(N)] with the last time t(N)=t0 t(k)<t(k-1)
t0=t(end);s=t0-t;
if(rows(s)>1) s=flipud(s);
else s=fliplr(s);
end
G=@(x,s) -F(x,t0-s);
X=lsode(G,x0,s);
X=flipud(X);

end

X=backivp(@(x,t) x,1,t=linspace(-1,0));
plot(t,X,t,exp(t),";exp(t);");pause(1);



X=backivp(@(x,t) t,1,t=linspace(-1,0));
%x(t)=x(0)-\int_t^0 s ds= 1+0.5*t^2
plot(t,X,t,1+0.5*t.*t);pause(1);


%graph of 1. s --> x(5,s) on [-1,3] while x(t,s) solves x'=x x(0)=s
%         2. s --> x(1,s) on [-2,2] while x(t,s) solves x''=-10sin(x) x(0)=0 x'(0)=s

s=linspace(-1,3);
xs=zeros(size(s));
for k=1:length(s),
  X=lsode(@(x,t) x,s(k),[0,5]);
  xs(k)=X(2);
end  
plot(s,xs);



s=linspace(-2,2);
xs=zeros(size(s));
for k=1:length(s),
  X=lsode(@(x,t) [x(2);-10*sin(x(1))],[0,s(k)],[0,1]);
  xs(k)=X(2,1);
end  
plot(s,xs);



