%
%  Risolve il problema BVP
%
%  y'' + p(x) y' + q(x) y = r(x)
%
%  con diferenze centrate al secondo ordine.
%  Il sistema risultante
% 
%  +                                            + +        +   +        +
%  |  alpha(1)  gamma(1)                        | | y(1)   |   | b(1)   |
%  |  beta(2)   alpha(2)  gamma(2)              | |        |   |        |
%  |              .......                       | |  ..    | = |        |
%  |                                            | |        |   |        |
%  |                       beta(N-1) alpha(N-1) | | y(N-1) |   | b(N-1) |
%  +                                            + +        +   +        +
%
%  E' risolto con le primitive di MATLAB
%
function [x,y]=BVPsolve(a,ya,b,yb,N)
  h = (b-a)/N ;
  %for i=1:N-1
  %  x(i) = a + i*h ;
  %end
  x    = [a+h:h:b-h] ;
  % riempimento del vettore alpha
  al   = -2 + (h^2).*q(x) ;
  % riempimento del vettore beta
  be   = 1 - h.*p(x)./2 ;
  % riempimento del vettore gamma
  ga   = 1 + h.*p(x)./2 ;
  % riempimento del vettore r
  bvec = h^2.*r(x) ;
  % condizioni al contorno
  bvec(1)   = bvec(1)   - be(1)   * ya ; 
  bvec(N-1) = bvec(N-1) - ga(N-1) * yb ; 
  %
  A = sparse(diag(al(1:N-1),0)) + ... 
      sparse(diag(be(2:N-1),-1)) + ...
      sparse(diag(ga(1:N-2),+1)) ;
  %A = diag(al(1:N-1),0)  + ...
  %	  diag(be(2:N-1),-1) + ...
  %	  diag(ga(1:N-2),+1) ;
  % soluzione del sistema
  y = A\bvec' ;
  spy(A)
end
