function exercise3 ( )

%*****************************************************************************80
%
%% exercise3 prints the solution of the tiny Poisson problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise3:\n' );
  fprintf ( 1, '  Print the solution of the tiny Poisson problem.\n' );

  a = 0.0;
  b = 3.0;
  c = 0.0;
  d = 2.0;
  nx = 4;
  hx = ( b - a ) / ( nx - 1 );
  x = linspace ( a, b, nx );
  ny = 3;
  hy = ( d - c ) / ( ny - 1 );
  y = linspace ( c, d, ny );

  A = zeros(nx*ny,nx*ny);  
  rhs = zeros(nx*ny,1);

  k = 0;
  for j = 1 : ny
    for i = 1 : nx
      k = k + 1;
      if ( j == 1 || j == ny || i == 1 || i == nx )
        A(k,k) = 1.0;
        rhs(k) = g ( x(i), y(j) );
      else
        CENTER = k;
        A(k,CENTER) = 2.0 / hx^2 + 2.0 / hy^2;
        WEST = k_from_ij ( i-1, j, nx );
        A(k,WEST) = -1.0 / hx^2;
        EAST = k_from_ij ( i+1, j, nx );
        A(k,EAST) = -1.0 / hx^2;
        SOUTH = k_from_ij ( i, j-1, nx );
        A(k,SOUTH) = -1.0 / hy^2;
        NORTH = k_from_ij ( i, j+1, nx );
        A(k,NORTH) = -1.0 / hy^2;
        rhs(k) = f ( x(i), y(j) );
      end
    end
  end
  
  u = A \ rhs;
%
%  Tabulate the solution.
%
  fprintf ( 1, ' k   i   j    x       y       u       exact\n' );
  k = 0;  
  for j = 1 : ny
    fprintf ( 1, '\n' );
    for i = 1 : nx
      k = k + 1;
      fprintf ( 1, '%2d  %2d  %2d  %6.2f  %6.2f  %6.2f  %6.2f\n', ...
                      k,   i,   j, x(i), y(j), u(k), g(x(i),y(j)) );
    end
  end

  return
end

function k = k_from_ij ( i, j, nx )

%*****************************************************************************80
%
  k = i + ( j - 1 ) * nx;
  return
end

function value = f ( x, y )

%*****************************************************************************80
%
  value = - 2.0 * y;
  return
end

function value = g ( x, y )

%*****************************************************************************80
%
  value = x.^2 .* y;
  return
end

