function [] = prob_3_6(X)
% 
% 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 Omega gamma;
Omega = 1;
gamma = 0.5; 

solinit = bvpinit(linspace(0,X),@guess);  
sol = bvp4c(@f,@bcs,solinit);
figure; 
p1 = plot( sol.x, sol.y(1,:), '-bo' ); hold on; 
legend( [ p1 ], 'y(x)' ); 
grid on; title( [ 'X = ', num2str(X) ] );


function [ v ]  = guess(x)
% GUESS - 
%
global Omega gamma;  
  
ept = exp(-Omega*x); 
v = [ (1/Omega)*(1-ept); ept; -Omega*ept; (Omega^2)*ept ]; 


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


function [ ode ] = f(x,y)
% F - 
% y(1)=f(z); y(2)=f'(z); y(3)=f''(z); y(4)=f'''(z)
%  
global Omega gamma;  

ode = [  
    y(2); 
    y(3);
    y(4);
    (-Omega+y(1))*y(4) + gamma*(y(2)*y(3) + Omega*(y(2))^2);
      ];