function prob_3_27
% 
% 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.
% 
%-----
  
% y(1) = y(x), y(2) = z(x), y(3) = y'(x), y(4) = z'(x)
global d
%options = [];
options = bvpset('Vectorized','on'); 
hold on
for i = 1:3
	D = 4*i;
    d = 1/D;
	if i == 1
        P = [-1; 0.5];  % Guess parameters P = [b; c]:
        solinit = bvpinit(linspace(d,D,5),[d; 0; 0; 0],P);
	else
		solinit = bvpinit(sol,[d,D]);
	end
	sol = bvp4c(@odes,@bcs,solinit,options);
	plot(sol.x,sol.y(1:2,:),sol.x(end),sol.y(1:2,end),'ro',...
         sol.x(1),sol.y(1:2,1),'ko');
    legend('y(x)','z(x)',0);
	axis([0 12 0 1]);
	drawnow
end
hold off

%===========================================================
function dydx = odes(x,y,P)
dydx = [ y(3,:); y(4,:); 2*(y(3,:) - y(2,:)) ./ (3*y(1,:)); -y(1,:).*y(4,:) ];
 
function res = bcs(ya,yb,P)
global d
res = zeros(6,1);
res(1:4) = ya - series(d,P);
res(5:6) = [ yb(3); (yb(4) - (- yb(1)*yb(2))) ];

function y = series(x,P)
b = P(1);
c = P(2);
yx  = x + b*x^(5/3) - c*x^2 - (5/7)*b^2*x^(7/3) ...
      + (7/6)*b*c*x^(8/3);
ypx = 1 + (5/3)*b*x^(2/3) - 2*c*x - (5/3)*b^2*x^(4/3)...
      + (28/9)*b*c*x^(5/3);
zx  = 1 + c*x - (1/6)*c*x^3 - (9/88)*b*c*x^(11/3);
zpx = c - (1/2)*c*x^2 - (3/8)*b*c*x^(8/3);
y = [ yx; zx; ypx; zpx];