%
% Solve the BVP of the form
%
%   y'' + p y' + q y = r
%   y(a) = ya, y(b) = yb
%
% with N subinterval and the shooting method
%
%   h    = (b-a)/N
%   x(k) = a + h*k,   
%   y(k) ~ y(a+h*k),
%
%  p : reference of a function of one argument p(x)
%  q : reference of a function of one argument q(x)
%  r : reference of a function of one argument r(x)
%
%  BCa : vector of two element [ a , ya ] with left BC
%  BCb : vector of two element [ b , yb ] with ight BC
%
%  N : number of subdivision of interval [a,b]
%
function [x,y] = BVPshoot( p, q, r, BCa, BCb, N )

  global pFunction qFunction rFunction ; % store problem in global variables

  pFunction = p ;
  qFunction = q ;
  rFunction = r ;

  % extract the boundary conditions
  a  = BCa(1) ;
  b  = BCb(1) ;
  ya = BCa(2) ;
  yb = BCb(2) ;
  % fill x
  h  = (b-a)/N ;
  x  = a + h*[0:N]' ;
  % first problem y(a) = ya, y'(a) = 0
  [x1,z] = ode45(@localOde1,x,[ya;0]) ;
  
  % second problem y(a) = 0, y'(a) = 1 and with r=0
  [x2,w] = ode45(@localOde2,x,[0;1]) ;
  
  % lambda for the linear combination solution
  lambda = (yb - z(end,1))/w(end,1) ; 
  
  % fill the y
  y = z(:,1)+lambda*w(:,1) ;
  
end

%
%  Transformation of the second oder ODE
%  
%   y''(x) + p(x) y'(x) + q(x) y(x) = r(x)
%  
%  To a first order ODE (system)
%
%   y'(x) = v(x)
%   v'(x) = r(x) - p(x) v(x) - q(x) y(x) 
%
%                    / y(x) \
%  The vector Z(x) = |      |
%                    \ v(x) /
%
function dZ = localOde1( x, Z )
  global pFunction qFunction rFunction ;

  dZ = [ Z(2) ; ...
         feval(rFunction,x) ...
         - feval(pFunction,x) *Z(2) ...
         - feval(qFunction,x) * Z(1) ] ;

end

%
%  Transformation of the second oder ODE
%  
%   y''(x) + p(x) y'(x) + q(x) y(x) = 0
%  
%  To a first order ODE (system)
%
%   y'(x) = v(x)
%   v'(x) = - p(x) v(x) - q(x) y(x) 
%
%                    / y(x) \
%  The vector Z(x) = |      |
%                    \ v(x) /
%
function dZ = localOde2( x, Z )
  global pFunction qFunction rFunction ;

  dZ = [ Z(2) ; ...
         - feval(pFunction,x) *Z(2) ...
         - feval(qFunction,x) * Z(1) ] ;

end
