function exercise2 ( )

%*****************************************************************************80
%
%% exercise2 sets up the linear system for the tiny discretized Poisson equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise2:\n' );
  fprintf ( 1, '  Set up the linear system for 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
  U = reshape ( u, nx, ny )
  [ X, Y ] = ndgrid ( x, y );
  G = g(X,Y)
  rms ( reshape ( U - G, nx*ny, 1 ) )

  return
end
function k = k_from_ij ( i, j, nx )

%*****************************************************************************80
%
%% k_from_ij gets the k coordinate from the (i,j) coordinates.
%
  k = i + ( j - 1 ) * nx;

  return
end
function value = g ( x, y )

%*****************************************************************************80
%
%% g returns the exact solution.
%
  value = x.^2 .* y;
  return
end
function value = f ( x, y )

%*****************************************************************************80
%
%% f returns the right hand side of the Poisson equation.
%
  value = - 2.0 * y;
  return
end
