function [u,x,y,error]=poissonfd(a,b,c,d,nx,ny,...
                         fun,bound,uex,varargin)
%POISSONFD approssimazione del problema di Poisson
%   in due dimensioni
%   [U,X,Y]=POISSONFD(A,B,C,D,NX,NY,FUN,BOUND)
%   risolve con lo schema alle differenze finite
%   a 5 punti il problema -LAPL(U) = FUN nel
%   rettangolo (A,B)X(C,D) con condizioni al bordo
%   di Dirichlet U(X,Y)=BOUND(X,Y) per ogni (X,Y)
%   sul bordo del rettangolo.
%   [U,X,Y,ERROR]=POISSONFD(A,B,C,D,NX,NY,...
%                 FUN,BOUND,UEX)
%   calcola anche l'errore nodale massimo ERROR
%   rispetto alla soluzione esatta UEX. FUN, BOUND
%   e UEX possono essere funzioni inline, anonymous
%   function o funzioni definite in M-file.
if nargin == 8
    uex = inline('0','x','y');
end
nx1 = nx+2; ny1=ny+2;
dim = nx1*ny1;
hx = (b-a)/(nx+1); hy = (d-c)/(ny+1);
    hx2 = hx^2;      hy2 = hy^2;
kii = 2/hx2+2/hy2; kix = -1/hx2;  kiy = -1/hy2;
K = speye(dim,dim);  rhs = zeros(dim,1);
y = c;
for m = 2:ny+1
 x = a; y = y + hy;
 for n = 2:nx+1
   i = n+(m-1)*nx1; x = x + hx;
   rhs(i) = feval(fun,x,y,varargin{:});
   K(i,i) = kii; K(i,i-1) = kix; K(i,i+1) = kix;
   K(i,i+nx1) = kiy;   K(i,i-nx1) = kiy;
 end
end
rhs1 = zeros(dim,1);
x = [a:hx:b]';
rhs1(1:nx1) = feval(bound,x,c,varargin{:});
rhs1(dim-nx-1:dim) = feval(bound,x,d,varargin{:});
y = [c:hy:d];
rhs1(1:nx1:dim-nx-1) = feval(bound,a,y,varargin{:});
rhs1(nx1:nx1:dim) = feval(bound,b,y,varargin{:});
rhs = rhs - K*rhs1;
nbound = [[1:nx1],[dim-nx-1:dim],[1:nx1:dim-nx-1],...
    [nx1:nx1:dim]];
ninternal = setdiff([1:dim],nbound);
K = K(ninternal,ninternal);
rhs = rhs(ninternal);
utemp = K\ rhs;
uh = rhs1; uh (ninternal) = utemp;
k = 1; y = c;
for j = 1:ny+1
    x = a;
    for i = 1:nx1
        u(i,j) = uh(k);         k = k + 1;
        ue(i,j) = feval(uex,x,y,varargin{:});
        x = x + hx;
    end
    y = y + hy;
end
x = [a:hx:b]; y = [c:hy:d];
if nargout == 4 & nargin == 8
        warning('Soluzione esatta non disponibile');
        error = [ ];
    else
        error = max(max(abs(u-ue)))/max(max(abs(ue)));
end
end