%systems_driver.m
% Solving systems of first order differential equations
% and multiple  plots       -    using MATLAB and ODE23.
% By Andy Ruina   (2/2/2000)

global  C1 C2  % This allows information to be passed to a function
               % without a direct pass.
 
C1 = 1;  C2 = -5;   % assign a value to the constants 

tzero  = 0;  % define the start value of t as 0
tfinal = 5;  %            final value of t as 5
             %  Hint:  make tfinal close to tzero when debugging,
			 %         then make it what you want when things work.
tspan = [tzero tfinal];  %this list will be given to ODE23

xzero = 2;  % The initial values for the integration
vzero = 3;
zzero  = [xzero vzero]; %this list will be given to ODE23

tol= 1e-4;  % a measure of the desired integration accuaracy.
            % a small number here will give a more accurate
			% solution but will take longer to run on the
			% computer.  1e-4 is a decent start.  But you can
			% use 1e-12 if you want a very accurate solution.
			% To get the computer to know this tolerance you
			% need the following line of code.
			
options = odeset('RelTol',tol,'AbsTol',tol);           

% 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', tspan, zzero, options);  % 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 160, 162, 163)

subplot(3,1,1)
 plot( t, z(:,1)) ; % 
 xlabel('time')    % this adds the text 'time' under the x axis
 ylabel('x')    
 title('Andy''s x vs t, v vs t, and state space v vs x:  2/2/2000') 
          
subplot(3,1,2)
 plot( t, z(:,2)) ; % 
 xlabel('time')    
 ylabel('v')    
 
subplot(3,1,3)
 plot(  z(:,1) , z(:,2) ) ;  %State space plot
 xlabel('x') 
 ylabel('v')