function prob_3_28 % % Written by: % -- % John L. Weatherwax 2005-05-07 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- global R color = [ 'k', 'r', 'b', 'k', 'r', 'b' ]; options = bvpset('FJacobian',@Jac,'BCJacobian',@BCJac,... 'Vectorized','on', 'Nmax', 50 ); %options = []; %options = bvpset( 'FJacobian', @Jac ); %options = bvpset( 'FJacobian', @Jac, 'BCJacobian', @BCJac ); %options = bvpset( 'Vectorized', 'on' ); R = 100; sol = bvpinit(linspace(0,1,5),ones(7,1),1); hold on for i = 1:5 sol = bvp4c(@ode,@bc,sol,options); fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters); fprintf('Number of mesh points = %f\n',length(sol.x)); plot(sol.x,sol.y(2,:),color(i)); axis([-0.1 1.1 0 1.7]); drawnow R = 10*R; end legend('R = 100','R = 1000','R = 10000',1); hold off %print -depsc ch3fig6 %================================================================ function dydx = ode(x,y,A); global R P = 0.7*R; dydx = [ y(2,:); y(3,:); R*(y(2,:).^2- y(1,:).*y(3,:)-A);... y(5,:); -R*y(1,:).*y(5,:)-1; y(7,:); -P*y(1,:).*y(7,:) ]; function [dFdy,dFdA] = Jac(x,y,A) global R dFdy = [ 0, 1, 0, 0, 0, 0, 0 0, 0, 1, 0, 0, 0, 0 -R*y(3), 2*R*y(2), -R*y(1), 0, 0, 0, 0 0, 0, 0, 0, 1, 0, 0 -R*y(5), 0, 0, 0, -R*y(1), 0, 0 0, 0, 0, 0, 0, 0, 1 -7/10*R*y(7), 0, 0, 0, 0, 0, -7/10*R*y(1) ]; dFdA = [ 0; 0; -R; 0; 0; 0; 0]; function res = bc(ya,yb,A) res = [ ya(1); ya(2); yb(1)-1; yb(2);... ya(4); yb(4); ya(6); yb(6)-1 ]; function [dBCdya,dBCdyb,dBCdA] = BCJac(ya,yb,A) dBCdya = [ 1, 0, 0, 0, 0, 0, 0 0, 1, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 1, 0, 0, 0 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 1, 0 0, 0, 0, 0, 0, 0, 0 ]; dBCdyb = [ 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0 1, 0, 0, 0, 0, 0, 0 0, 1, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 1, 0, 0, 0 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 1, 0 ]; dBCdA = zeros(8,1);