%
%  Solves Poisson's equation in a square with Dirichlet boundary
%  conditions:
%
%     u_xx + u_yy = f(x,y) ,  0 < x,y < 1
%
%     u(x,0) = g(x,0),  u(x,1) = g(x,1),  u(0,y) = g(0,y),  u(1,y) = g(1,y).
%
%  Uses a centered finite difference method.
%  I have set f(x,y) = x^2 + y^2 and g(x,y) = 1.

f = inline('x^2 + y^2','x','y');
g = inline('1','x','y');

%  Set up grid.

n = input('Enter number of subintervals in each direction: ');
h = 1/n;
N = (n-1)^2;

%  Initialize block tridiagonal finite difference matrix A and
%  right-hand side vector b.

A = sparse(zeros(N,N));  % This stores A as a sparse matrix.  Only the
                         % nonzero entries are stored.
b = zeros(N,1);

%  Loop over grid points in y direction.

for j=1:n-1,
  yj = j*h;

%    Loop over grid points in x direction.

  for i=1:n-1,
    xi = i*h;
    k = i + (j-1)*(n-1);    % k is the index of the equation corresponding
                            % to point (xi,yj)
    A(k,k) = -4/h^2;
    if i > 1, A(k,k-1) = 1/h^2; end;       % Coupling to point on left.
    if i < n-1, A(k,k+1) = 1/h^2; end;     % Coupling to point on right.
    if j > 1, A(k,k-(n-1)) = 1/h^2; end;   % Coupling to point below.
    if j < n-1, A(k,k+(n-1)) = 1/h^2; end; % Coupling to point above.

    b(k) = f(xi,yj);
    if i==1, b(k) = b(k) - g(0,yj)/h^2; end;   % Bndy point on left.
    if i==n-1, b(k) = b(k) - g(1,yj)/h^2; end; % Bndy point on right.
    if j==1, b(k) = b(k) - g(xi,0)/h^2; end;   % Bndy point below.
    if j==n-1, b(k) = b(k) - g(xi,1)/h^2; end; % Bndy point above.
  end;
end;

%  Solve finite difference equations.

u = A\b;

%  Plot results.

ugrid = reshape(u,n-1,n-1);  % This makes the Nx1 vector u into an n-1 x n-1
                             % array.
surf([h:h:1-h]',[h:h:1-h]',ugrid)  %This plots u(x,y) as a function of x and y.
title('Computed solution')