function [] = prob_3_3()
% PROB_3_3 - 
% 
% 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.
% 
%-----
  
d = 0.001; 
epsilon = 0.1; 
k = 0.1; 

solinit = bvpinit(linspace(d,1,10),@guess,0);  
sol = bvp4c(@f,@bcs,solinit,[],epsilon,k);
p = sol.parameters;
xsml = linspace(0,d,10);
p0=plot( xsml, p + (10/3)*(p/(0.1+p))*xsml.^2, '-og' ); hold on; 
p1=plot( sol.x, sol.y(1,:), '-bo' ); 
p2=plot( sol.x, sol.y(2,:), '-ro' ); grid on; 
legend( [p1,p2], 'y1', 'y2','Location','NorthWest' ); 
fprintf( 'p = %g.\n', sol.parameters );

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

v = [ y; yp ];


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


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

ode = [  y(2); 
	 y(1)/(epsilon*(y(1)+k)) - 2*y(2)/x ];