%
% Solve the BVP of the form
%
%   y'' + p y' + q y = r
%   y(a) = ya, y(b)=yb
%
% with N subinterval.
%
%   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 , y(a) ] with left BC
%  BCb : vector of two element [ b , y(b) ] with ight BC
%
%  N : number of subdivision of interval [a,b]
%
function [x,y] = BVPsolver( p, q, r, BCa, BCb, N )
  % 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]' ;
  % fill the 3 diagonal of the linear system and the RHS
  alpha = zeros( N-1, 1 ) ;
  beta  = zeros( N-1, 1 ) ;
  gamma = zeros( N-1, 1 ) ;
  omega = zeros( N-1, 1 ) ;
  
  alpha = -2/h^2 + feval( q, x(2:end-1) ) ;
  beta  = 1/h^2 - feval( p, x(2:end-1) ) ./(2*h) ;
  gamma = 1/h^2 + feval( p, x(2:end-1) ) ./(2*h) ;
  omega = feval( r, x(2:end-1) ) ;
  omega(1)   = omega(1)   - beta(1)*ya ;
  omega(N-1) = omega(N-1) - gamma(N-1)*yb ;

  % this is equivalent to the following loop
  % for k=1:N-1
  %   alpha(k) = -2/h^2 + feval( q, x(k+1) ) ;
  %   beta(k)  = 1/h^2 - feval( p, x(k+1) ) /(2*h) ;
  %   gamma(k) = 1/h^2 + feval( q, x(k+1) ) /(2*h) ;
  %   omega(k) = feval( r, x(2:end-1) ) ;
  % end
  % omega(1)   = omega(1)   - beta(1)*ya ;
  % omega(N-1) = omega(N-1) - gamma(N-1)*yb ;
  
  % fill the tridiagonal system
  A = zeros( N-1, N-1 ) ;
  A = diag(alpha,0) + diag(gamma(1:N-2),1) + diag( beta(2:N-1),-1) ;

  % this is equivalent to the following loop
  % for k=1:N-1
  %   A(k,k) = alpha(k) ;
  % end  
  % for k=1:N-2
  %   A(k+1,k) = beta(k+1) ; % fill the lower diagonal
  %   A(k,k+1) = gamma(k)  ; % fill the upper diagonal
  % end
  
  % fill the RHS
  b = zeros( N-1, 1 ) ;
  b = omega ;

  % solve the linear system
  sol = A\b ;
  
  % fill the y
  y = [ ya; sol; yb] ;

end