function [] = prob_3_1(lambda0)
% PROB_3_1 - 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-04-13
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

%solinit = bvpinit(linspace(0,1,10),@guess,10);
%solinit = bvpinit(linspace(0,1,10),@guess,100);
%solinit = bvpinit(linspace(0,1,10),@guess,500);
%solinit = bvpinit(linspace(0,1,10),@guess,1000);
%solinit = bvpinit(linspace(0,1,10),@guess,5000);
solinit = bvpinit(linspace(0,1,10),@guess,lambda0);  
sol = bvp4c(@f,@bcs,solinit);
p1=plot( sol.x, sol.y(1,:), '-bo' ); hold on; 
p2=plot( sol.x, sol.y(2,:), '-ro' ); grid on; 
legend( [p1,p2], 'y', 'v','Location','NorthWest' ); 
fprintf( 'lambda = %g.\n', sol.parameters );

function [res] = bcs(ya,yb,lambda)
% BCS - 
%   
res = [ ya(1); yb(1)-1; yb(2) ];


function [ v ]  = guess(x)
% GUESS - 
%   
y = x.*sin(2.5*pi*x);
yp = sin(2.5*pi*x) + 2.5*pi*x.*cos(2.5*pi*x);
v = [ y; (1-x.^7).*yp ];


function [ ode ] = f(x,y,lambda)
% F - 
% y(1)=y; y(2)=v
% 

% Incorporate the singular behaviour explicitly:
if( x > .9995 ) 
  ode = [ lambda/7;
	  -lambda*(x^7)*y(1) ];
else
  ode = [ y(2)/(1-x^7);
	  -lambda*(x^7)*y(1) ];
end