function [ u, it_num ] = poisson_1d ( n, a, b, ua, ub, force )

%*****************************************************************************80
%                                                    
%% poisson_1d() solves a 1D PDE, using the Gauss-Seidel method.
%
%  Discussion:
%
%    This routine solves a 1D boundary value problem of the form
%
%      - U''(X) = F(X) for A < X < B,
%
%    with boundary conditions U(A) = UA, U(B) = UB.
%
%    The Gauss-Seidel method is used. 
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 October 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Hager,
%    Applied Numerical Linear Algebra,
%    Prentice-Hall, 1988,
%    ISBN13: 978-0130412942,
%    LC: QA184.H33.
%
%  Input:
%
%    integer N, the number of intervals.
%
%    real A, B, the left and right endpoints of the region.
%
%    real UA, UB, the left and right boundary values.
%
%    function value = FORCE ( x ), the name of the function 
%    which evaluates the right hand side.
%
%  Output:
%
%    integer IT_NUM, the number of iterations.
%
%    real U(N+1,1), the computed solution.
%

%
%  Initialization.
%
  tol = 0.0001;
%
%  Set the nodes.
%
  x = ( linspace ( a, b, n + 1 ) )';
%
%  Set the right hand side.
%
  r = zeros ( n + 1, 1 );

  r(1) = ua;
  r(2:n) = force ( x(2:n) ) / n / n;
  r(n+1) = ub;
%
%  Initialize the solution.
%
  u = zeros ( n + 1, 1 );
  u(1) = ua;
  u(n+1) = ub;
%
%  Gauss-Seidel iteration.
%
  it_num = 0;

  while ( true )

    it_num = it_num + 1;

    [ u, d1 ] = gauss_seidel ( n + 1, r, u );

    if ( d1 <= tol )
      break
    end

  end

  return
end
