subroutine gauss_seidel ( n, r, u, dif_l1 )

c*********************************************************************72
c                                                    
cc gauss_seidel() carries out one step of a Gauss-Seidel iteration.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 November 2011
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    William Hager,
c    Applied Numerical Linear Algebra,
c    Prentice-Hall, 1988,
c    ISBN13: 978-0130412942,
c    LC: QA184.H33.
c
c  Input:
c
c    integer N, the number of unknowns.
c
c    double precision R(N), the right hand side.
c
c    double precision U(N), the estimated solution.
c
c  Output:
c
c    double precision U(N), the updated estimated solution.
c
c    double precision DIF_L1, the L1 norm of the difference between the
c    input and output solution estimates.
c
      implicit none

      integer n

      double precision dif_l1
      integer i
      double precision r(n)
      double precision u(n)
      double precision u_old

      dif_l1 = 0.0D+00
c
c  Setting U(1) = R(1), U(N) = R(N) seems right, but it makes the
c  code fail.
c
c     u(1) = r(1)
      do i = 2, n - 1
        u_old = u(i)
        u(i) = 0.5D+00 * ( u(i-1) + u(i+1) + r(i) )
        dif_l1 = dif_l1 + abs ( u(i) - u_old )
      end do
c     u(n) = r(n)

      return
      end
      subroutine poisson_1d ( n, a, b, ua, ub, force, it_num, u )

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

      integer n

      double precision a
      double precision b
      double precision d1
      external force
      double precision force
      double precision h
      integer i
      integer it_num
      double precision r(n+1)
      double precision tol
      double precision u(n+1)
      double precision ua
      double precision ub
      double precision x(n+1)
c
c  Initialization.
c
      tol = 0.0001D+00
c
c  Set the nodes.
c
      do i = 1, n + 1
        x(i) = ( ( n + 1 - i ) * a + ( i - 1 ) * b ) / n
      end do
c
c  Set the right hand side.
c
      r(1) = ua
      h = ( b - a ) / dble ( n )
      do i = 2, n
        r(i) = h**2 * force ( x(i) )
      end do
      r(n+1) = ub
c
c  Initialize the solution.
c
      u(1) = ua
      do i = 1, n
        u(i) = 0.0D+00
      end do
      u(n+1) = ub
c
c  Gauss-Seidel iteration.
c
      it_num = 0

10    continue

      it_num = it_num + 1

      call gauss_seidel ( n + 1, r, u, d1 )

      if ( tol < d1 ) then
        go to 10
      end if

      return
      end