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