function prob_3_30_pt2
% 
% 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.
% 
%-----
  
mu = 0.02; 
lambda = 0.0279;
eta = 0.01; 
beta0V = [ 10, 100, 500, 1000 1575 ]; 

options = bvpset( 'Vectorized', 'on' ); 
for bi=1:length(beta0V),
  beta0=beta0V(bi);
  if( bi==1 )
    sol = bvpinit(linspace(0,1,10),[1 1 1]);
  else
    sol = bvpinit(sol,[0,1]); 
  end
  sol = bvp4c(@odes,@bcs,sol,options,mu,lambda,eta,beta0);
  figure; plot(sol.x,sol.y(1,:),...
	       sol.x,100*sol.y(2,:),...
	       sol.x,100*sol.y(3,:));
  title( [ 'Beta0=',num2str(beta0) ] );
end




%=================================================
function dydx = odes(x,y,mu,lambda,eta,beta0)
  
beta = beta0*(1+cos(2*pi*x));

dydx = [ mu - beta .* y(1,:) .* y(3,:);
         beta .* y(1,:) .* y(3,:) - y(2,:)/lambda;
         y(2,:)/lambda - y(3,:)/eta;
	 ];

function res = bcs(ya,yb,mu,lambda,eta,beta0)
res = [ ya(1)-yb(1); 
	ya(2)-yb(2);
	ya(3)-yb(3); 
	];