% % 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 k2x0    = 0;x_end = 10;y0    = [ 1; 0; 0];k1    = 1;k2    = 10;rhs_fctn   = 'ChemReac';h = 0.2;% Exact solutionmx    = ceil((x_end-x0)/ h);x     = (x0:h:x0+mx*h);[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,1)plot(x,y_exp(1,:),'g',x,y_exp(2,:),'r',x,y_exp(3,:),'b');legend('y_1', 'y_2', 'y_3') xlabel(' x ')title([' Explicit Euler, h = ', num2str(h) ])[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,3)plot(x,y_imp(1,:),'g',x,y_imp(2,:),'r',x,y_imp(3,:),'b');legend('y_1', 'y_2', 'y_3') xlabel(' x ')title([' Implicit Euler, h = ', num2str(h) ])% Compute exact solutionA = [ -k1   0   0;       k1  -k2  0;       0    k2  0 ];[V,D] = eig(A);y_ex = zeros(3,length(x));for j = 1:length(x)    y_ex(:,j)  = V*(exp(diag(D)*x(j)).*(V\y0));end% Plot the error in the computed solutionsubplot(2,2,2)plot(x,abs(y_exp(1,:)-y_ex(1,:)),'g'); hold onplot(x,abs(y_exp(2,:)-y_ex(2,:)),'r'); hold onplot(x,abs(y_exp(3,:)-y_ex(3,:)),'b'); hold offlegend('y_1', 'y_2', 'y_3') xlabel(' x ')title([' Error, Explicit Euler, h = ', num2str(h) ])% Plot the error in the computed solutionsubplot(2,2,4)plot(x,abs(y_imp(1,:)-y_ex(1,:)),'g'); hold onplot(x,abs(y_imp(2,:)-y_ex(2,:)),'r'); hold onplot(x,abs(y_imp(3,:)-y_ex(3,:)),'b'); hold offlegend('y_1', 'y_2', 'y_3') xlabel(' x ')title([' Error, Implicit Euler, h = ', num2str(h) ])returndisp(' Hit key to continue .....')pause%  Solution of the %                         %    y_1'(x) =  y_2(x)%               %    y_2'(x) =  4*y_1(x) + 3*y_2(x) - 8*exp(x) * cos(2x)%                %  with initial conditions%%    y_1(0) =  2 + 10/13             %    y_2(0) = -3 + 14/13%%  This system is equivalent to the second order ODE%%    z''(x) - 3 z'(x) - 4 z(x) = - 8 exp(x) * cos(2x)%%  The exact solution of the system is given by%    y_1(x) =  exp(-x) + exp(4*x) %               + 10/13 * exp(x) * cos(2x) + 2/13 * exp(x) * sin(2x)%               %    y_2(x) = -exp(-x) + 4 * exp(4*x) %               + 14/13 * exp(x) * cos(2x) - 18/13 * exp(x) * sin(2x)%clf;x0    = 0;x_end = 1;y0    = [ 2 + 10/13; 3 + 14/13];rhs_fctn   = 'Ex1';h = 0.02;% Exact solutionmx    = ceil(x_end-x0)/ h;x     = (x0:h:x0+mx*h);y_ex  = [ exp(-x) + exp(4*x) + (10/13)*exp(x).*cos(2*x) ...              + (2/13)*exp(x).*sin(2*x);         -exp(-x) + 4*exp(4*x) + (14/13)*exp(x).*cos(2*x) ...              - (18/13).*exp(x).*sin(2*x) ];[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,1)plot(x,y_exp(1,:),'g'); hold onplot(x,y_exp(2,:),'r'); hold offlegend('y_1', 'y_2') xlabel(' x ')title([' Explicit Euler, h = ', num2str(h) ])% Plot the error in the computed solutionsubplot(2,2,2)plot(x,abs(y_exp(1,:)-y_ex(1,:)),'g'); hold onplot(x,abs(y_exp(2,:)-y_ex(2,:)),'r'); hold offlegend('y_1', 'y_2') xlabel(' x ')title([' Error, Explicit Euler, h = ', num2str(h) ])[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,3)plot(x,y_imp(1,:),'g'); hold onplot(x,y_imp(2,:),'r'); hold offlegend('y_1', 'y_2') xlabel(' x ')title([' Implicit Euler, h = ', num2str(h) ])% Plot the error in the computed solutionsubplot(2,2,4)plot(x,abs(y_imp(1,:)-y_ex(1,:)),'g'); hold onplot(x,abs(y_imp(2,:)-y_ex(2,:)),'r'); hold offlegend('y_1', 'y_2') xlabel(' x ')title([' Error, Implicit Euler, h = ', num2str(h) ])disp(' Hit key to continue .....')pause%  Predator Prey model 1%%%    y_1'(x) = alpha*y_1(x) - beta*y_1(x)*y_2(x)%    y_2'(x) = gamma*y_2(x)  + delta*y_1(x)*y_2(x)%%%clf;x0    = 0;x_end = 20;y0    = [ 80; 30];rhs_fctn   = 'PredPrey';h = 0.05;[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,1)plot(x,y_exp(1,:),'g'); hold onplot(x,y_exp(2,:),'r'); hold offlegend('Population 1', 'Population 2') xlabel(' Time ')title([' Explicit Euler, h = ', num2str(h) ])[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,2)plot(x,y_imp(1,:),'g'); hold onplot(x,y_imp(2,:),'r'); hold offlegend('Population 1', 'Population 2') xlabel(' Time ')title([' Implicit Euler, h = ',  num2str(h) ])% Plot the difference in computed solutionssubplot(2,2,3)plot(x,abs(y_exp(1,:)-y_imp(1,:)),'g'); hold onplot(x,abs(y_exp(2,:)-y_imp(2,:)),'r'); hold offlegend('Population 1', 'Population 2') xlabel(' Time ')title([' Difference in Computed Solutions, h = ',  num2str(h) ])disp(' Hit key to continue .....')pause%  Predator Prey model 2%%                                         y_1(x) * y_2(x)%    y_1'(x) = 1.2 * y_1(x) - y_1^2(x) -  ---------------%                                          y_1(x) + 0.2%%               1.5 * y_1(x) * y_2(x)%    y_2'(x) =  ---------------------  -  y_2(x)%                   y_1(x) + 0.2%%%clf;x0    = 0;x_end = 30;y0    = [ 1; 0.75];rhs_fctn   = 'PredPrey2';h = 0.001;[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,1)plot(x,y_exp(1,:),'g'); hold onplot(x,y_exp(2,:),'r'); hold offlegend('Population 1', 'Population 2') xlabel(' Time ')title([' Explicit Euler, h = ', num2str(h) ])[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );% Plot the computed solutionsubplot(2,2,2)plot(x,y_imp(1,:),'g'); hold onplot(x,y_imp(2,:),'r'); hold offlegend('Population 1', 'Population 2') xlabel(' Time ')title([' Implicit Euler, h = ',  num2str(h) ])% Plot the difference in computed solutionssubplot(2,2,3)plot(x,abs(y_exp(1,:)-y_imp(1,:)),'g'); hold onplot(x,abs(y_exp(2,:)-y_imp(2,:)),'r'); hold offlegend('Population 1', 'Population 2') xlabel(' Time ')title([' Difference in Computed Solutions, h = ',  num2str(h) ])