fprintf ( 1, '\n' );
fprintf ( 1, 'poisson_test\n' );
fprintf ( 1, '  Test poisson(), a simple poisson equation solver\n' );
fprintf ( 1, '  for a rectangular region.\n' );

nx = 21;
ny = 41;

xmin = 0.0;
xmax = 1.0;
ymin = 1.0;
ymax = 3.0;

[ U, X, Y ] = poisson ( nx, ny, xmin, xmax, ymin, ymax, @f, @g );
%
%  Assuming g(x,y) is the exact solution, compute the RMS error.
%  You could also print the computed and exact solutions, for comparison.
%
err = 0.0;
for i = 1 : ny
  for j = 1 : nx
    err = err + ( U(i,j) - g(X(i,j),Y(i,j)) )^2;
  end
end
err = sqrt ( err / nx / ny );

fprintf ( 1, '\n' );
fprintf ( 1, '  Using a %d by %d grid\n', nx, ny );
fprintf ( 1, '  The RMS error was %g\n', err );
%
%  Plot the computed and exact solutions.
%
figure ( 1 )
surf ( X, Y, U );
title ( 'Computed solution' );
print ( '-dpng', 'poisson_computed_surface.png' );

figure ( 2 )
U2 = g(X,Y);
surf ( X, Y, U2 );
title ( 'Exact solution' );
print ( '-dpng', 'poisson_exact_surface.png' );

figure ( 3 )
contourf ( X, Y, U );
title ( 'Computed solution' );
print ( '-dpng', 'poisson_computed_contours.png' );

function value = f ( x, y )
%
%% f evaluates the right hand side of the Poisson equation.
%
  value = pi * pi * ( x .* x + y .* y ) .* sin ( pi * x .* y );
  return
end
function value = g ( x, y )
%
%% g evaluates the boundary conditions.
%
%  g is only needed for the boundary conditions.
%  But if it is actually a formula for the exact solution everywhere,
%  it can also be used to compute the computational error.
%
  value = sin ( pi * x .* y );
  return
end
