function prob_3_30_pt1
% 
% 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; 
beta0 = 1575;

options = bvpset( 'Vectorized', 'on' ); 
sol = bvpinit(linspace(0,1,10),(10^(-3))*[1 1 1 1 1 1]);
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,:));

%=================================================
function dydx = odes(x,y,mu,lambda,eta,beta0)
  
sz = size(x); 
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;
	 zeros(sz);
	 zeros(sz);
	 zeros(sz); 
	 ];

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