function [xh,uh,xi,f]=bvpmod(a,b,N,mu,eta,sigma,bvpfun,...
                    ua,ub,varargin)
%BVP Risolve un problema ai limiti
%  [XH,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.
%  XH e UH contengono, rispettivamente, i nodi
%  e la soluzione numerica, inclusi i valori al bordo.
h = (b-a)/(N+1);
xh = (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 = xh(2:end-1);
%f =bvpfun(xi,varargin{:});
f =bvpfun;
f(1) =  f(1)+ua*(hm+hd);
f(end) = f(end)+ub*(hm-hd);
uh = A\f;
uh=[ua; uh; ub];