% Simple scalar ODE-IVP y'= y+exp(t)*cos(t), y(0)= 0.
% Define the ODE function.
fun = @(t,y)y+exp(t)*cos(t);
% Set the options to print out the run statistics. 
options = odeset('Stats','on');
% Call ode45 with tspan = [0 pi], y(0) = 0, and the options. 
[T,Y] = ode45(fun,[0 pi],0,options);
% Plot the results.
plot(T,Y,'-b',T,Y,'r.');figure(1);


%%
% Simple oscillator with drag, released from rest with unit
% displacement at time t = 0. 
% The second-order ODE is x'' + Dx + x = 0, where x is displacement and D
% is a drag coefficient. The initial conditions are x(0) = 1, x'(0) = 0. 
% Rewritten as an equivalent first order system to feed to ode45, this is 
% 
%     y_1' = y_2, y_2' = -y_1 - Dy_2,    
%     y_1(0) = 1, y_2(0) = 0, 
%
% where y_1 = x, y_2 = x'. In matrix-vector form, which ode45 wants this is 
%
%                     (1)
%     y' = Ay, y(0) = (0)
%
%           (y_1)         (0   1)
% where y = (y_2) and A = (-1 -D). (I'm doing my best to represent
% 2-vectors and 2x2 matrices here.) 
%
% Set D as desired. Try D = 0 (no drag) and D = .1, 1, 2. 
D=0; 
% Set A for this D.  
A=[0 1;-1 -D];
% Define the ODE function.
fun = @(t,y)A*y;
% Set the options to print out the run statistics. 
options = odeset('Stats','on');
% Call ode45 with tspan = [0 10*pi], y(0) = [1;0], and the options. 
[T,Y] = ode45(fun,[0 10*pi],[1;0],options);
% Uncomment for solution-vs-time plot:
plot(T,Y(:,1),T,Y(:,2));figure(1);
% Uncomment for phase plane plot:
%plot(Y(:,1),Y(:,2));figure(1);axis square;


%%
% Predator-prey model in which y_1 and y_2 are the populations of prey and
% predators (think rabbits and foxes). 
% Set alpha.
alpha = .01;
% Define the ODE function. 
fun = @(t,y)[2*y(1)-alpha*y(1)*y(2);-y(2)+alpha*y(1)*y(2)];
% Set the options to print out the run statistics. 
options = odeset('Stats','on','RelTol',1.e-6,'AbsTol',1.e-9);
% Call ode45 with tspan = [0 20], y(0) = [300;100], and the options. 
[T,Y] = ode45(fun,[0 20],[300;100],options);
% Uncomment for solution-vs-time plot:
plot(T,Y(:,1),T,Y(:,2));figure(1);
% Uncomment for phase plane plot:
%plot(Y(:,1),Y(:,2));figure(1);axis square;
