function main1
options= odeset ‘RelTol’,1e-6,‘Stats’,‘on’ ;
%initial conditions
T=100;
Xo = 1; 0; T+30 ;
tspan = 0,10 ;
tic
t,X = ode45 @TestFunction,tspan,Xo,options ;
toc
figure
hold on
plot t, X :,1
plot t, X(:,2
plot t, X :,3 ,‘.’
legend ‘x1’,‘x2’,‘x3’
ylabel ‘x’
xlabel ‘t’
return
function dx_dt = TestFunction t,x
k1=2;k2=0.1;k3=10 -3 ;k4=2.66*10 -3 ;K=4.44*10 -2 ;T=100;
dx_dt(1) =-k1*x(1)+k2*x(2)2-2*k3*x(1)2;
dx_dt(2) =2*k1*x(1)-2*(k2+k4)*x(2)2;
dx_dt(3)=K*(T-x(3));
dx_dt = dx_dt’;
return
Algorithm 1: Scilab/Matlab program to find the numerical solution of (2) and (3).