function [] = prob_3_2()
% 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.
% 
%-----
  
global d;  
d = 0.1; 

solinit = bvpinit(linspace(d,16,10),@guess,1);  
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], 'y1', 'y2','Location','NorthWest' ); 
fprintf( 'p = %g.\n', sol.parameters );

function [res] = bcs(ya,yb,p)
% BCS - 
%   
global d;   
  
res = [ ya(1)-(1/10)-p*sqrt(d)-(1/10^3)*d; 
	ya(2)-(1/2)*p*sqrt(10)-(1/10^3); 
	yb(1)-(1/6) ];


function [ v ]  = guess(x)
% GUESS - 
%
% $$$ y = zeros(1,length(x));
% $$$ yp = y;

% WWX: Alternativly one could use
% (as suggested in the text):
p = 1; 
y = (1/10) + p*sqrt(x) + (1/10^3)*x;
yp = y; 

v = [ y; yp ];


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

% Incorporate the singular behaviour explicitly:
ode = [  y(2); 
	(y(1)^3 - y(2))/(2*x); ];