% Solving systems of first order differential equations
% and multiple plots - using MATLAB and ODE23.
% By Andy Ruina (9/30/96)
% This is a script file (see Pratap pgs 29-38).
% It is VERY similar to the sample used for a single first order ODE.
global C1 C2 % This allows information to be passed to a function
% without a direct pass.
C1 = 2; C2 = .5; % assign a value to the constants
tzero = 2; % define the start value of t as 2
tfinal = 20; % final value of t as 20
% Hint: make tfinal close to tzero when debugging,
% then make it what you want when things work.
zzero = [1.5 4]; % The solution z is a now a list of numbers.
% This sets the initial condition as z1(tzero) = 1.5,
% and z2(tzero) = 4.
tol = 1e-4; % tol is an optional argument of ODE23 that sets
% the accuracy of the solution.
% Hint: make tol big when debugging, then small for final
% accurate answer.
% the actual solution is done in the line that follows.
% The result of this command is two arrays: t and z, where
% t contains a list of times (chosen by the computer but with
% the first element tzero and the last tfinal) and z an array
% of values of z at those times. Each row of z corresponds to
% a given time. Each element in the row the value of z1 or z2.
% For example z(8,2) is the value of z2 at the eighth time step.
% type >help ode23 at the MATLAB prompt,
% or see Pratap pgs 82-83.
[t z] = ode23 ('systems_rhs', tzero, tfinal, zzero, tol ); % the key line
%To get multiple plots on one page, each with its own label and
%axis, use sublot (type >help subplot or see Pratap pages 114, 116, 117)
subplot(3,1,1)
plot( t, z(:,1)) ; %
xlabel('time') % this adds the text 'time' under the x axis
ylabel('z(1)')
title('Andys z1 vs t, z2 vs t, and state space z2 vs z1: 9/30/96')
subplot(3,1,2)
plot( t, z(:,2)) ; %
xlabel('time')
ylabel('z(2)')
subplot(3,1,3)
plot( z(:,1) , z(:,2) ) ; %State space plot
xlabel('z1')
ylabel('z2')
% Note this is not quite a 'phase plane' plot because the system is
% 'non-autonomous', the right hand side depends on time. In a 2D
% phase plane plot with no explicit time dependence on the right hand
% side, solutions cannot cross (why?).