function [] = prob_3_21
% 
% 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, [0.5]);
sol = bvp4c(@f,@bcs,solinit,[]);

% Plot t v.s. y1(t):figure;
t0 = sol.parameters;
plot( (1-t0)*sol.x+t0, sol.y(1,:), '-bo' ); grid on;
fprintf( 'reduced drag coefficient = %f\n', sol.y(3,1) );


function [v] = guess(x)

v = [ x; 1; 1 ];

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

function [ ode ] = f(x,y,P)
% F - 
% y(1)=y(x); y(2)=y'(x); y(3)=t^2 + \int_t^1 \frac{2 \tau}{ y'(\tau)^2 + 1 }
% 

t0 = P(1);

ode = [      (1-t0)*y(2);
	     ( (1-t0)/(t0+(1-t0)*x) ) * ( y(2)*(1+(y(2)^2))/(3*(y(2)^2)-1) );
	     2*(1-t0)*(t0+(1-t0)*x)*(y(2)^2) / ( (y(2)^2) + 1 );
      ];