%  Uses the five-point formula to solve Poisson's equation
%  on a square.

%  Inline functions define the source function (f), the boundary
%  values (ux0,ux1,u0y,u1y), and the true solution (u).

f = inline('x^2 + y^2');
ux0 = inline('0');
ux1 = inline('.5*x^2');
u0y = inline('sin(pi*y)');
u1y = inline('exp(pi)*sin(pi*y) + .5*y^2');
u = inline('exp(pi*x)*sin(pi*y) + .5*(x^2)*(y^2)');

%  Ask user for problem size.  (Use same no. of points in each direction.)

n = input(' Enter number of subintervals in each direction: ');
h = 1/n;
N = (n-1)^2;
A = sparse(N,N);   % Store A as a sparse matrix.
F = zeros(N,1);

%  Set up matrix A and right-hand side vector F.

A = -4*sparse(eye(N,N));
for j=1:n-1,                % Loop over grid rows
  for i=1:n-1,                % Loop over points in each row

    k = (j-1)*(n-1)+i;          % Index of this point
    if j > 1, 
      A(k,k-(n-1)) = 1;                    % coupling to pt below
    end;
    if j < n-1,
      A(k,k+(n-1)) = 1;                    % coupling to pt above
    end;
    if i > 1, A(k,k-1) = 1; end;           % coupling to pt to left
    if i < n-1, A(k,k+1) = 1; end;         % coupling to pt to right

    xi = i*h;  yj = j*h;
    F(k) = f(xi,yj);                         % right-hand side vector
    if j==1, 
      F(k) = F(k) - (1/h^2)*ux0(xi);         % bdry pt below
    end;
    if j==n-1,
      F(k) = F(k) - (1/h^2)*ux1(xi);         % bdry pt above
    end;
    if i==1,
      F(k) = F(k) - (1/h^2)*u0y(yj);         % bdry pt to left
    end;
    if i==n-1,
      F(k) = F(k) - (1/h^2)*u1y(yj);         % bdry pt to right
    end;

  end;
end;

%  Remember to multiply A by 1/h^2.

A = (1/h^2)*A;

%  Solve linear system.

uapprox = A\F;

%  Compare to true solution.

utrue = zeros(N,1);
for j=1:n-1,
  for i=1:n-1,
    k = (j-1)*(n-1)+i;
    xi = i*h;  yj = j*h;
    utrue(k) = u(xi,yj);
  end;
end;
err2 = h*norm(utrue-uapprox),
errinf = norm(utrue-uapprox,'inf')

%  Plot solution and error over grid.
ugrid = reshape(utrue,n-1,n-1);
figure(1)
mesh([h:h:(n-1)*h]',[h:h:(n-1)*h]',ugrid')
title('True Solution')
errgrid = reshape(utrue-uapprox,n-1,n-1);
figure(2)
mesh([h:h:(n-1)*h]',[h:h:(n-1)*h]',errgrid')
title('Error')