% 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?).