function sol = prob_3_29
% 
% 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.
% 
%-----
  
%options = [];  
options = bvpset('Vectorized','on'); 

epsilonV = 1 ./ [ 10 20 30 50 75 100 ];  
for epsi=1:length(epsilonV),
  epsilon=epsilonV(epsi);
  if epsi==1
    solinit = bvpinit(linspace(-pi/2,pi/2,20),0.5,1);
  else
    solinit = bvpinit(sol,[-pi/2,pi/2]);
  end
  sol = bvp4c(@ode,@bcs,solinit,options,epsilon);

  fprintf('lambda = %g.\n',sol.parameters);
  figure; plot(sol.x,sol.y(1,:));
  axis([-pi/2 pi/2 0 1.1]); drawnow;
end

%=================================================
function dydx = ode(x,y,lambda,epsilon)
%dydx = (sin(x)^2 - lambda*(sin(x)^4)/y)/epsilon;
dydx = (sin(x).^2 - lambda*(sin(x).^4)./y(1,:))/epsilon;

function res = bcs(ya,yb,lambda,epsilon)
res = [ ya-1; yb-1 ];