%
% Genero la matrice dei coefficienti e il temine noto del 
% problema ellittico
%
%  u(x,y)    + u(x,y)      = f(x,y)     su (0,1) x (0,1)
%        xx          yy
%
%  u(x,y) = g(x,y)  su bordo (0,1) x (0,1)
%
%  N = numero i suddivisioni dell'intervallo [0,1]
%  f = funzione termine sorgente
%  g = funzione che implementa le condizioni al contorno
%
%  A = matrice dei coefficienti
%  b = termine noto
%
%  la soluzione � posta nelle metrici X, Y, e Z che contengono
%  i valori nodali nella mesh compreso i termini al bordo
%
%  X = X(i,j) coodinata x del nodo (i,j)
%  Y = Y(i,j) coodinata y del nodo (i,j)
%  Z = Z(i,j) coodinata z della soluzione nel nodo (i,j)
%
%
%         +                 +
%         |        1        |
%    1    |                 |
%  -----  |  1    -4     1  |  u_ij = f_ij
%   h^2   |                 |
%         |        1        |
%         +                 +
%
function [A,b,X,Y,Z] = elliptic(N,f,g)
  %
  h = 1/N ; % dimensione della mesh
  
  % numerazione delle equazioni
  EQ = zeros(N+1,N+1) ;
  
  % numero le equazioni (solo punti interni)
  EQ(2:end-1,2:end-1) = [1:N-1]'* ones(1,N-1) + ...
                        ones(N-1,1)*[0:N-2]*(N-1) ;
  neq = (N-1)^2 ;
  % questo corrisponde al ciclo
  % neq = 0 ;
  % for i=2:N
  %   for j=2:N
  %     neq     = neq + 1 ;
  %     EQ(i,j) = neq ;
  %   end
  % end

  % ordinamento RED-BLACK
  % neq = 0 ;
  % for i=2:N
  %   for j=2:N
  %     if mod(i+j,2) == 0
  %       neq     = neq + 1 ;
  %       EQ(i,j) = neq ;
  %     end
  %   end
  % end
  % for i=2:N
  %   for j=2:N
  %     if mod(i+j,2) == 1
  %       neq     = neq + 1 ;
  %       EQ(i,j) = neq ;
  %     end
  %   end
  % end

  % coefficienti dello stencil
  h2 = h^2 ;

  % costruzione della matrice dei coefficienti
  A = sparse( neq, neq ) ;
  b = zeros ( neq, 1   ) ;
  for i=2:N
    x = (i-1)*h ;
    for j=2:N
      y = (j-1)*h ;

      % numero equazioni vicine
      ic = EQ(i,j) ;
      il = EQ(i-1,j) ; % left
      ir = EQ(i+1,j) ; % right
      it = EQ(i,j+1) ; % top
      ib = EQ(i,j-1) ; % bottom

      % termine noto
      b(ic) = h2*feval( f, x, y ) ;

      % coefficienti dello stencil
      A(ic,ic) = -4 ;
      
      if il > 0
        % Punto interno aggiungo alla matrice il pezzo di stencil
        A(ic,il) = 1 ;
      else
        % Punto di bordo aggiungo al termine noto la BC
        b(ic) = b(ic) - feval( g, x-h, y ) ;      
      end

      if ir > 0
        A(ic,ir) = 1 ;
      else
        b(ic) = b(ic) - feval( g, x+h, y ) ;      
      end
      
      if it > 0
        A(ic,it) = 1 ;
      else
        b(ic) = b(ic) - feval( g, x, y+h ) ;      
      end

      if ib > 0
        A(ic,ib) = 1 ;
      else
        b(ic) = b(ic) - feval( g, x, y-h ) ;   
      end
    end
  end
  
  % calcolo la soluzione
  SOL = A\b ;

  % riempimento delle matrici in uscita, inizializzazione
  Z     = zeros(N+1,N+1) ;
  [X,Y] = meshgrid(linspace(0,1,N+1),linspace(0,1,N+1)) ;
  
  % Cerco gli indici dei nodi interni, quelli numerati
  IDX    = find( EQ > 0 ) ;

  % copio la soluzione tenendo conto dell'ordine delle equazioni
  Z(IDX) = SOL(EQ(IDX)) ;

  % Cerco gli indici dei nodi di bordo
  IDX    = find( EQ == 0 ) ;

  % copio i valori al contorno sulla soluzione
  Z(IDX) = feval( g, X(IDX), Y(IDX) ) ;

end
