function [u,x,y,error]=poissondf(a,b,c,d,nx,ny,f,g,uex)
%POISSONDF Il metodo a 5 punti per il problema di Poisson.
%   [U,X,Y]=POISSONDF(A,B,C,D,NX,NY,F,G) calcola la soluzione
%   approssimata U del problema di Poisson con termine
%   forzante F e dato di Dirichlet G sul rettangolo
%   [A,B]X[C,D] su una griglia di (NX+1) x (NY+1) nodi
%   le cui coordinate sono restituite nelle matrici X e Y.
% 
%   [U,X,Y,ERROR]=POISSONDF(A,B,C,D,NX,NY,F,G,UEX) restituisce
%   in ERROR l'errore calcolato in norma infinito discreto
%   rispetto alla soluzione esatta UEX.
%
%   Si veda anche VALLIM.

%   Saleri Fausto 24-01-01

if nargin == 8 
  uex = '0 + 0.*x'; 
end

hx = (b-a)/nx; 
hy = (d-c)/ny; 
nx1 = nx+1;
hx2 = hx^2; 
hy2 = hy^2; 
kii = 2/hx2+2/hy2; 
kix = -1/hx2; 
kiy = -1/hy2;
dim = (nx+1)*(ny+1);
K = speye(dim,dim); 
rhs = zeros(dim,1);
y = c;
for m = 2:ny
   x = a; y = y + hy;
   for n = 2:nx
      i = n+(m-1)*(nx+1);
      x = x + hx; 
      rhs(i) = eval(f);  
      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);
y = c; 
x = [a:hx:b]; 
rhs1(1:nx1) = eval(g);
y = d; 
rhs1(dim-nx:dim) = eval(g);
x = a; 
y = [c:hy:d]; 
rhs1(1:nx1:dim-nx) = eval(g);
x = b; 
rhs1(nx1:nx1:dim) = eval(g);
rhs = rhs - K*rhs1;
nbordo = [[1:nx1],[dim-nx:dim],[1:nx1:dim-nx],[nx1:nx1:dim]];
ninterni = setdiff([1:dim],nbordo);
K = K(ninterni,ninterni); 
rhs = rhs(ninterni); 
utemp = K\rhs; 
uh = rhs1; 
uh (ninterni) = 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) = eval(uex); 
     x = x + hx;
  end
  y = y + hy;
end
x = [a:hx:b]; 
y = [c:hy:d];
if nargout == 4 
   if nargin == 8
      error = 'Manca la soluzione esatta';
   else
      error = max(max(abs(u-ue)));
   end 
end

return
