کد:
% Start fresh ----------------------------------------------------------
clear all
% The following line passes Alpha and Beta into a function
global Alpha Beta
% Program header -------------------------------------------------------
disp('Simulate a non-isothermal 1st-order reactor: ')
disp(' dc/dt = 1 - c - xk*c ... material balance')
disp(' dT/dt = 1 - T - alpha*xk*c ... enthalpy balance')
disp(' ')
% Read/assign model parameters -----------------------------------------
Alpha = -5.;
Beta = 2.;
% Initial values (row vector) ------------------------------------------
y0(1) = 1. ;
y0(2) = 0.1;
% Integration limits ---------------------------------------------------
t0 = 0.; % Start time
tf = 1.; % Finish time
% Call a routine to solve ODE ------------------------------------------
% [t y] = ode45('odefcn', t0, tf, y0); %MATLAB v4
[t y] = ode45('odefcn', [t0:0.01:tf], y0); %MATLAB v5
% Save the result in a file
temp = [t y]; save ode.dat temp -ascii
% save ode.dat [t y] -ascii %This does not work!
% Plot the result
plot(t, y)
% plot(t, 10*y(:,1), t, y(:,2)) %Plot with different scales
% plot(y(:,1), y(:,2)) % plot y1 versus y2 (phase trajectory)
% The following command saves the plot in a file named ODE.MET
% meta ode %obsolute in v.4 use "print"
% To convert ODE.MET into a postscript file name OUTPUT.PS, issue from DOS
% DOS>matprint ode
% In the basement of Eng. Class Bldg, send OUTPUT.PS to a printer from Windows