function ch3ex5
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' ];
color = [ 'k','k', 'k', 'k', 'k' ];
wave_speed = [5, 10, 15, 20, 25];
hold on
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));
    if i == 1
       plot(sol.x,sol.y(1,:),'-k');
    end 
    if i == 2
       plot(sol.x,sol.y(1,:),':k');
    end   
    if i == 3
       plot(sol.x,sol.y(1,:),'--k');
    end    
    if i == 4
       plot(sol.x,sol.y(1,:),'-.k');
    end    
    if i == 5
       plot(sol.x,sol.y(1,:),'-k');
    end   
    axis([-250 250 0 1]);
    drawnow
end
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 ];