function poisson_2d_test ( )

%*****************************************************************************80
%
%% poisson_2d_test() tests poisson_2d().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 October 2024
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../poisson_2d' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'poisson_2d_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test poisson_2d().\n' );
%
%  Problem 1.
%
  nx = 11;
  ny = 11;
  poisson_2d ( nx, ny, @rhs1, @u_exact1, @uxxyy_exact1 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'poisson_2d_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../poisson_2d' )

  return
end
function f = rhs1 ( nx, ny )

%*****************************************************************************80
%
%% rhs1() initializes the right hand side "vector" for problem 1.
%
%  Discussion:
%
%    It is convenient for us to set up RHS as a 2D array.  However, each
%    entry of RHS is really the right hand side of a linear system of the
%    form
%
%      A * U = F
%
%    In cases where U(I,J) is a boundary value, then the equation is simply
%
%      U(I,J) = F(i,j)
%
%    and F(I,J) holds the boundary data.
%
%    Otherwise, the equation has the form
%
%      (1/DX^2) * ( U(I+1,J)+U(I-1,J)+U(I,J-1)+U(I,J+1)-4*U(I,J) ) = F(I,J)
%
%    where DX is the spacing and F(I,J) is the value at X(I), Y(J) of
%
%      pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 September 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, NY, the X and Y grid dimensions.
%
%  Output:
%
%    real F(NX,NY), the right hand side data.
%
  f = zeros ( nx, ny );

  for j = 1 : ny
    y = ( j - 1 ) / ( ny - 1 );
    for i = 1 : nx
      x = ( i - 1 ) / ( nx - 1 );
%
%  The "boundary" entries of F store the boundary values of the solution.
%
      if ( i == 1 || i == nx || j == 1 || j == ny )
        f(i,j) = u_exact1 ( x, y );
%
%  The "interior" entries of F store the right hand sides of the Poisson equation.
%
      else
        f(i,j) = - uxxyy_exact1 ( x, y );
      end
    end
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function u = u_exact1 ( x, y )

%*****************************************************************************80
%
%% uexact1() evaluates the exact solution for problem 1.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, Y, the coordinates of a point.
%
%  Output:
%
%    real U, the value of the exact solution at (X,Y).
%
  u = sin ( pi * x .* y );

  return
end
function uxxyy = uxxyy_exact1 ( x, y )

%*****************************************************************************80
%
%% uxxyy_exact1() evaluates exact ( d/dx d/dx + d/dy d/dy ) for problem 1.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, Y, the coordinates of a point.
%
%  Output:
%
%    real UXXYY, the value of ( d/dx d/dx + d/dy d/dy ) of the
%    exact solution at (X,Y).
%
  uxxyy = - pi * pi * ( x .* x + y .* y ) .* sin ( pi * x .* y );

  return
end
 
