function [] = prob_3_19
% 
% 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.
% 
%-----

epsilon=1;
alpha = 0.5;

solinit = bvpinit(linspace(0,1,10),@guess, [1]);  
sol = bvp4c(@f,@bcs,solinit,[],alpha,epsilon);

figure; 
plot( sol.x, sol.y(1,:), '-bo' ); hold on; grid on;
disp( 'Values requested' ); 
disp( [ 'omega = ', num2str(sol.parameters) ] ); 
deval( sol, [ 0.5, 1 ] )

function v = guess(x)
% GUESS - 
%   
v = [ cos(pi*x); -pi*sin(pi*x); x; 1 ];


function [res] = bcs(ya,yb,P,alpha,epsilon)
% BCS - 
%      

res = [ ya(2);
        ya(1)-1;
        ya(3);
        yb(2);
	yb(4)-(1/alpha^2)*(1-(1-alpha^2)*yb(3));
      ];


function [ ode ] = f(x,y,P,alpha,epsilon)
% F - 
% y(1)=y(x); y(2)=y'(x); y(3)=\int_0^x ...; y(4)=(1/alpha^2)[...
% 
  
omega = P(1);  

tmp = 1/sqrt( 1 + (epsilon^2)*(y(1)^2) );
  
ode = [  
    y(2);
    -(omega^2)*( ((1-alpha^2)/y(4)) * tmp + alpha^2 )*y(1);
    tmp; 
    0; 
      ];