function exercise1 ( )

%*****************************************************************************80
%
%% exercise1 solves the steady heat equation in a square, with gaussian solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise1:\n' );
  fprintf ( 1, '  Solve the steady heat equation in a square.\n' );

  nx = 11;
  a = -1.0;
  b = +1.0;
  x = linspace ( a, b, nx );
  dx = ( b - a ) / ( nx - 1 );

  ny = 11;
  c = -1.0;
  d = +1.0;
  y = linspace ( c, d, ny );
  dy = ( d - c ) / ( ny - 1 );

  [ X, Y ] = ndgrid ( x, y );
  A = zeros(nx*ny,nx*ny);  
  rhs = zeros(nx*ny,1);

  for j = 1 : ny
    for i = 1 : nx

      [ C, N, S, E, W ] = cnsew ( i, j, nx );

      if ( j == 1 | j == ny | i == 1 | i == nx )
        A(C,C) = 1.0;
        rhs(C) = g ( x(i), y(j) );
      else
        A(C,C) =  2.0 / dx^2 + 2.0 / dy^2;
        A(C,W) = -1.0 / dx^2;
        A(C,E) = -1.0 / dx^2;
        A(C,S) =              -1.0 / dy^2;
        A(C,N) =              -1.0 / dy^2;
        rhs(C) = f ( x(i), y(j) );
      end
    end
  end
  
  u = A \ rhs;
  U = reshape ( u, nx, ny );
%
%  Since we know the exact solution, we can print the RMS error.
%
  V = g ( X, Y );
  E = rms ( reshape ( U - V, nx*ny, 1 ) );
  fprintf ( 1, '  RMS error = %g\n', E );
%
%  Plot the solution.
%

  surf ( X, Y, U );
  title ( 'Solution of steady heat equation with source term.' );
  colorbar ( );

  print ( '-dpng', 'exercise1.png' );

  return
end

function [ C, N, S, E, W ] = cnsew ( i, j, nx )

%*****************************************************************************80
%
  C = i     + ( j - 1 ) * nx;
  N = i     + ( j     ) * nx;
  S = i     + ( j - 2 ) * nx;
  E = i + 1 + ( j - 1 ) * nx;
  W = i - 1 + ( j - 1 ) * nx;
  return
end

function value = g ( x, y )

%*****************************************************************************80
%
  value = exp ( - x.^2 - y.^2 );
  return
end

function value = f ( x, y )

%*****************************************************************************80
%
  value = ( 4.0 - 4.0 * x.^2 - 4.0 * y.^2 ) * exp ( - x.^2 - y.^2 );
  return
end
