%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           Solving odes explicitly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% ODEs can sometimes be solved explicitly
   %%% with the command 'dsolve'. (see 'help dsolve' for details)
   %%% Example: y'(t) = t*y(t)
y = dsolve('Dy = t*y','t')
   % Answer: y = C1*exp(t^2/2)


   %%% With initial value: y'(t) = t*y(t), y(0) = 2
y = dsolve('Dy=t*y','y(0)=2','t')
   % Answer: y = 2*exp(t^2/2)

   %%% The returned object 'y' is actually a symbolic object
   %%% and we cannot plot it directly. For this reason we need
   %%% to convert it into a valid function. First we vectorize()
   %%% it so that it also accepts vectors as input.
   %%% Then we plot it at 100 points of the interval [0,1] 
t = linspace(0,1,100);
z = eval(vectorize(y));
hold off
plot(t,z)
input('Press return for next figure')


   %%% Next example: y'(t)=y(t)*(100-y(t))
y = dsolve('Dy=y*(100-y)')
   % Answer: y=
   %         0
   %         100
   %         -100/(exp(C2-100*t)-1)

   %%% Higher order example: y''(t)+8y'(t)+2y(t)=cos(t), y(0)=0, y'(0)=1
y = dsolve('D2y+8*Dy+2*y=cos(t)','y(0)=0,Dy(0)=1')
t = linspace(0,1,100);
z = eval(vectorize(y));
hold off
plot(t,z)
input('Press return for next figure')

   %%% Multi-dimensional example: x'(t)=x(t)+2y(t),  x(0)=1
   %%%                            y'(t)=4x(t)-y(t),  y(0)=0
[x,y] = dsolve('Dx=x+2*y','Dy=x-4*y','x(0)=1,y(0)=0')
t = linspace(0,1,100);
u = eval(vectorize(x));
v = eval(vectorize(y));
hold off
plot(t,u)
hold on
p = plot(t,v)
set(p,'Color','red')
input('Press return for next figure')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           Solving odes numerically
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% ODEs can be solved numerically
   %%% with the command 'ode45'. (see 'help ode45' for details)
   %%% Example: y'(t) = t*y(t), y(0)=2
   %%% First define the right-hand side
f = inline('t*y')
   %%% Syntax:  ode45(function,domain,initial condition)
[t,y] = ode45(f,[0,3],2);
   %%% (Note: Matlab chooses the discretization step size by itself)
hold off
plot(t,y)

   %%% For systems of odes, it is convenient to use M-files.
   %%% Let the content of the file 'firstode.m' be as follows:
%%% function yprime = firstode(x,y);
%%% yprime = t*y;
   %%% Then the above code works if the inline function is replaced
   %%% by a pointer @ to indicate the M-file:
[t,y] = ode45(@firstode,[0,3],2);
hold off
plot(t,y)
input('Press return for next figure')

   %%% Next Example: Lotka-Volterra Predator-Prey
   %%% prey:     x'(t) = 3*x(t)-x(t)*y(t)      x(0)=10
   %%% predator: y'(t) = 0.1*x(t)*y(t)-2*y(t)  y(0)=1
%%% function xprime = Lotka(t,y);
%%% xprime = [3*x(1)-x(1)*x(2);0.1*x(1)*x(2)-2*x(2)];
   %%% Note that x is stored as x(1) and y is stored as x(2)
   %%% Note that the semicolons indicate that xprime is a column vector
[t,x] = ode45(@Lotka,[0,20],[10,1]);
hold off
plot(t,x(:,1))
hold on
p=plot(t,x(:,2));
set(p,'Color','red')
   %% Next plot the trajectory:
input('Press return for next figure')
hold off
plot(x(:,1),x(:,2))
