function prob_3_25
% 
% 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 c alpha beta infty
options = []; 
options = bvpset(options,'Vectorized','on');
options = bvpset(options,'FJacobian',@odeJac);
options = bvpset(options,'BCJacobian',@bcJac);
color = [ 'r', 'b', 'g', 'm', 'k' ];
wave_speed = [5, 10, 15, 20, 25];
hold on
tic; 
for i = 1:5
	c = wave_speed(i);
	alpha = (-c + sqrt(c^2 + 4))/2;
    beta  = (-c + sqrt(c^2 - 4))/2;
    infty = 10*c;
%    if i == 1
        solinit = bvpinit(linspace(-infty,infty,20),@guess);
%    else
%        solinit = bvpinit(sol,[-infty,infty]);
%    end
    sol = bvp4c(@ode,@bc,solinit,options);
    plot(sol.x,sol.y(1,:),color(i));
    axis([-250 250 0 1]);
    drawnow
end
toc
legend('c = 5', 'c = 10', 'c = 15', 'c = 20', 'c = 25',3);
hold off
%print -depsc ch3fig5
%======================================================
function v = guess(z)
global c alpha beta infty
if z > 0
    v = [exp(beta*z); beta*exp(beta*z)];
else
    v = [(1 - exp(alpha*z)); -alpha*exp(alpha*z)];
end

function dydz = ode(z,y)
global c alpha beta infty
dydz =  [ y(2,:); -(c*y(2,:) + y(1,:).*(1 - y(1,:))) ];

function dFdy = odeJac(z,y)
global c alpha beta infty
dFdy = [       0,          1
         (-1 + 2*y(1)),  - c ];

function res = bc(ya,yb)
global c alpha beta infty
res = [ ya(2)/(ya(1) - 1) - alpha
        yb(1)/exp(beta*infty) - 1 ];

% deriveFbc.m derives these partial derivatives.
function [dBCdya, dBCdyb] = bcJac(ya,yb)
global c alpha beta infty
dBCdya = [ -ya(2)/(ya(1) - 1)^2, 1/(ya(1) - 1)
                   0                   0       ];

dBCdyb = [         0          0
           1/exp(beta*infty), 0 ];