% % ode_drv.m: simple driver for ode solvers%%%  Solution of the chemical reaction system%                         %    y_1'(x) = -k1*y_1(x) %    y_2'(x) =  k1*y_1(x) - k2*y_2(x)%    y_3'(x) =              k2*y_2(x)%%clf;% declare reaction rates as global variablesglobal k1 k2y0    = [ 1; 0; 0];k1    = 1;k2    = 10;rhs_fctn   = 'ChemReac23';reltol = 1.e-3;abstol = 1.e-3;opts   = odeset('reltol',reltol,'abstol', abstol);[t, y] = ode23tx( rhs_fctn, [0,10], y0, opts );% Plot the computed solutionfigure(1)plot(t,y(:,1),'g+',t,y(:,2),'r+',t,y(:,3),'b+');legend('y_1', 'y_2', 'y_3') xlabel(' t ')title(['ode23tx (abstol = ', num2str(abstol), 'reltol = ', num2str(reltol),')'])% Compute exact solutionA = [ -k1   0   0;       k1  -k2  0;       0    k2  0 ];[V,D] = eig(A);y_ex = zeros(length(t),3);for j = 1:length(t)    y_ex(j,:)  = (V*(exp(diag(D)*t(j)).*(V\y0)))';end% Plot the error in the computed solutionfigure(2)semilogy(t,abs(y(:,1)-y_ex(:,1)),'g'); hold onsemilogy(t,abs(y(:,2)-y_ex(:,2)),'r'); hold onsemilogy(t,abs(y(:,3)-y_ex(:,3)),'b'); hold offlegend('y_1', 'y_2', 'y_3') xlabel(' t ')title(['ode23tx (abstol = ', num2str(abstol), 'reltol = ', num2str(reltol),')'])