function [x,uh]=bvp(a,b,N,mu,eta,sigma,bvpfun,...
                    ua,ub,varargin)
%BVP Risolve un problema ai limiti
%   [X,UH]=BVP(A,B,N,MU,ETA,SIGMA,BVPFUN,UA,UB) risolve
%   con il metodo delle differenze finite centrate in
%   N nodi equispaziati interni ad (A,B) il problema
%     -MU*D(DU/DX)/DX+ETA*DU/DX+SIGMA*U=BVPFUN
%   sull'intervallo (A,B) con condizioni al bordo
%   U(A)=UA e U(B)=UB.
%   BVPFUN puo' essere una funzione inline, una
%   anonymous function o una funzione definita in
%   un M-file.
%   X e UH contengono, rispettivamente, i nodi
%   e la soluzione numerica, inclusi i valori al bordo.
h = (b-a)/(N+1);
x = (linspace(a,b,N+2))';
hm = mu/h^2;  hd = eta/(2*h);
e =ones(N,1);
A = spdiags([-hm*e-hd (2*hm+sigma)*e -hm*e+hd],...
    -1:1, N, N);
xi = x(2:end-1);
f =feval(bvpfun,xi,varargin{:});
f(1) =  f(1)+ua*(hm+hd);
f(end) = f(end)+ub*(hm-hd);
uh = A\f;
uh=[ua; uh; ub];
