program main

!*****************************************************************************80
!
!! MAIN is the main program for FEM1D_NONLINEAR.
!
!  Discussion:
!
!    FEM1D_NONLINLEAR solves a nonlinear one dimensional boundary value problem.
!
!    The differential equation has the form:
!
!      -d/dx (p(x) du/dx) + q(x)*u +u*u' =  f(x)
!
!    The finite-element method uses piecewise linear basis functions.
!
!    Here U is an unknown scalar function of X defined on the
!    interval [XL,XR], and P, Q and F are given functions of X.
!
!    The values of U or U' at XL and XR are also specified.
!
!    Sample problem #1:
!
!    u(x)  = x,
!    p(x)  = 1,
!    q(x)  = 0,
!    f(x)  = x,
!    u(0)  = 0,
!    u'(1) = 1.
!    The code should solve this problem exactly.
!
!    Sample problem #2:
!
!    u(x)  = 2*(1-cos(0.5*pi*x))/pi,
!    p(x)  = 1,
!    q(x)  = 0,
!    f(x)  = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0)  = 0,
!    u'(1) = 1.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    29 April 2007
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    real ( kind = rk ) ADIAG(NU).
!    ADIAG(I) is the "diagonal" coefficient of the I-th
!    equation in the linear system.  That is, ADIAG(I) is
!    the coefficient of the I-th unknown in the I-th equation.
!
!    real ( kind = rk ) ALEFT(NU).
!    ALEFT(I) is the "left hand" coefficient of the I-th
!    equation in the linear system.  That is, ALEFT(I) is the
!    coefficient of the (I-1)-th unknown in the I-th equation.
!    There is no value in ALEFT(1), since the first equation
!    does not refer to a "0-th" unknown.
!
!    real ( kind = rk ) ARITE(NU).
!    ARITE(I) is the "right hand" coefficient of the I-th
!    equation in the linear system.  ARITE(I) is the coefficient
!    of the (I+1)-th unknown in the I-th equation.  There is
!    no value in ARITE(NU) because the NU-th equation does not
!    refer to an "NU+1"-th unknown.
!
!    real ( kind = rk ) F(NU).
!    ASSEMBLE stores into F the right hand side of the linear
!    equations.
!    SOLVE replaces those values of F by the solution of the
!    linear equations.
!
!    real ( kind = rk ) FOLD(NU).
!    FOLD contains the value of F from the previous iteration,
!    and is used in ASSEMBLE to add correction terms to the
!    matrix and right hand side.
!
!    real ( kind = rk ) H(N), the length of the subintervals.
!
!    integer IBC.
!    IBC declares what the boundary conditions are.
!    1, at the left endpoint, U has the value UL,
!       at the right endpoint, U' has the value UR.
!    2, at the left endpoint, U' has the value UL,
!       at the right endpoint, U has the value UR.
!    3, at the left endpoint, U has the value UL,
!       and at the right endpoint, U has the value UR.
!    4, at the left endpoint, U' has the value UL,
!       at the right endpoint U' has the value UR.
!
!    integer IMAX.
!    The number of Newton iterations to carry out.
!
!    integer INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    integer NL.
!    The number of basis functions used in a single
!    subinterval.  (NL-1) is the degree of the polynomials
!    used.  For this code, NL is fixed at 2, meaning that
!    piecewise linear functions are used as the basis.
!
!    integer NODE(NL,N).
!    For each subinterval I:
!    NODE(1,I) is the number of the left node, and
!    NODE(2,I) is the number of the right node.
!
!    integer NPRINT.
!    The number of points at which the computed solution
!    should be printed out when compared to the exact solution.
!
!    integer NQUAD.
!    The number of quadrature points used in a subinterval.
!    This code uses NQUAD = 1.
!
!    integer N.
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    integer NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    real ( kind = rk ) UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    real ( kind = rk ) UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    real ( kind = rk ) XL.
!    XL is the left endpoint of the interval over which the
!    differential equation is being solved.
!
!    real ( kind = rk ) XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
!    real ( kind = rk ) XQUAD(N)
!    XQUAD(I) is the location of the single quadrature point
!    in interval I.
!
!    real ( kind = rk ) XR.
!    XR is the right endpoint of the interval over which the
!    differential equation is being solved.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: n = 10
  integer, parameter :: nl = 2

  real ( kind = rk ) adiag(n+1)
  real ( kind = rk ) aleft(n+1)
  real ( kind = rk ) arite(n+1)
  real ( kind = rk ) f(n+1)
  real ( kind = rk ) h(n)
  integer i
  integer ibc
  integer imax
  integer indx(0:n)
  integer node(nl,n)
  integer nprint
  integer nquad
  integer nu
  integer problem
  real ( kind = rk ) u(n+1)
  real ( kind = rk ) ul
  real ( kind = rk ) ur
  real ( kind = rk ) xl
  real ( kind = rk ) xn(0:n)
  real ( kind = rk ) xquad(n)
  real ( kind = rk ) xr

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEM1D_NONLINEAR'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Solve a nonlinear boundary value problem:'
  write ( *, '(a)' ) '    -d/dx (p(x) du/dx) + q(x)*u + u*u'' =  f(x)'
  write ( *, '(a)' ) '  on an interval [xl,xr], with the values of'
  write ( *, '(a)' ) '  u or u'' specified at xl and xr.'
!
!  Do two problems
!
  do problem = 1, 2
!
!  Initialize variables that define the problem.
!
    call init ( ibc, imax, nprint, nquad, problem, ul, ur, xl, xr )
!
!  Compute the quantities that describe the geometry of the problem.
!
    call geometry ( h, ibc, indx, nl, node, n, nu, xl, xn, xquad, xr )

    u(1:nu) = 0.0
!
!  Begin the Newton iteration.
!
    do i = 1, imax
!
!  Assemble the linear system.
!
      if ( i <= 3 ) then
        call assemble_picard ( adiag, aleft, arite, f, h, indx, n, nl, &
          node, nquad, nu, problem, u, ul, ur, xn, xquad )
      else
        call assemble_newton ( adiag, aleft, arite, f, h, indx, n, nl, &
          node, nquad, nu, problem, u, ul, ur, xn, xquad )
      end if
!
!  Print out the linear system, just once.
!
      if ( i == 1 ) then
        call prsys ( adiag, aleft, arite, f, nu )
      end if
!
!  Solve the linear system.
!
      call solve ( adiag, aleft, arite, f, nu, u )
!
!  Print the current solution.
!
      call output ( u, ibc, indx, n, nu, ul, ur, xn )

    end do
!
!  Compare the solution to the exact solution.
!
    call compare ( u, indx, n, nl, node, nprint, nu, problem, ul, &
      ur, xl, xn, xr )

  end do
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEM1D_NONLINEAR:'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine assemble_newton ( adiag, aleft, arite, f, h, indx, n, nl, &
  node, nquad, nu, problem, u, ul, ur, xn, xquad )

!*****************************************************************************80
!
!! ASSEMBLE_NEWTON assembles the Newton linear system.
!
!  Discussion:
!
!    The linear system being solved here is for the Newton correction
!    to an approximate solution of a nonlinear system.
!
!    Thus, we suppose that we have a nonlinear function F(X),
!    and an approximate solution X0.  If we can suppose there is an
!    exact solution X* that is "nearby", and in fact close enough
!    that Taylor's theorem gives us a useful estimate, then we
!    may write:
!
!      F(X*) = F(X0) + F'(X0) * ( X* - X0 ) + Order ( X* - X0 )^2
!
!    and by rearranging, we get the Newton system (which is only
!    approximately correct):
!
!      F'(X0) * ( X* - X0 ) = - F(X0)
!
!    We solve this system and add the solution to X0 to get a
!    new approximate solution that, we hope, is much improved.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    29 April 2007
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = rk ) ADIAG(NU).
!    ADIAG(I) is the "diagonal" coefficient of the I-th
!    equation in the linear system.  That is, ADIAG(I) is
!    the coefficient of the I-th unknown in the I-th equation.
!
!    Output, real ( kind = rk ) ALEFT(NU).
!    ALEFT(I) is the "left hand" coefficient of the I-th
!    equation in the linear system.  That is, ALEFT(I) is the
!    coefficient of the (I-1)-th unknown in the I-th equation.
!    There is no value in ALEFT(1), since the first equation
!    does not refer to a "0-th" unknown.
!
!    Output, real ( kind = rk ) ARITE(NU).
!    ARITE(I) is the "right hand" coefficient of the I-th
!    equation in the linear system.  ARITE(I) is the coefficient
!    of the (I+1)-th unknown in the I-th equation.  There is
!    no value in ARITE(NU) because the NU-th equation does not
!    refer to an "NU+1"-th unknown.
!
!    Output, real ( kind = rk ) F(NU), the right hand side of the linear
!    equations.
!
!    Input, real ( kind = rk ) H(N), the length of the subintervals.
!
!    Input, integer INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    Input, integer N.
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    Input, integer NL.
!    The number of basis functions used in a single
!    subinterval.  (NL-1) is the degree of the polynomials
!    used.  For this code, NL is fixed at 2, meaning that
!    piecewise linear functions are used as the basis.
!
!    Input, integer NODE(NL,N).
!    For each subinterval I:
!    NODE(1,I) is the number of the left node, and
!    NODE(2,I) is the number of the right node.
!
!    Input, integer NQUAD.
!    The number of quadrature points used in a subinterval.
!    This code uses NQUAD = 1.
!
!    Input, integer NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Input, real ( kind = rk ) U(NU), the solution value
!    from the previous iteration,
!
!    Input, real ( kind = rk ) UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    Input, real ( kind = rk ) UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    Input, real ( kind = rk ) XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
!    Input, real ( kind = rk ) XQUAD(N)
!    XQUAD(I) is the location of the single quadrature point
!    in interval I.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer nu

  real ( kind = rk ) aij
  real ( kind = rk ) adiag(nu)
  real ( kind = rk ) aleft(nu)
  real ( kind = rk ) arite(nu)
  real ( kind = rk ) f(nu)
  real ( kind = rk ) ff
  real ( kind = rk ) h(n)
  real ( kind = rk ) he
  integer ie
  integer ig
  integer il
  integer indx(0:n)
  integer iu
  integer iul
  integer iur
  integer iq
  integer jg
  integer jl
  integer jr
  integer ju
  integer nl
  integer node(nl,n)
  integer nquad
  real ( kind = rk ) phii
  real ( kind = rk ) phiix
  real ( kind = rk ) phij
  real ( kind = rk ) phijx
  real ( kind = rk ) pp
  integer problem
  real ( kind = rk ) qq
  real ( kind = rk ) total
  real ( kind = rk ) u(nu)
  real ( kind = rk ) ul
  real ( kind = rk ) uold
  real ( kind = rk ) uoldx
  real ( kind = rk ) ur
  real ( kind = rk ) x
  real ( kind = rk ) xleft
  real ( kind = rk ) xn(0:n)
  real ( kind = rk ) xqe
  real ( kind = rk ) xquad(n)
  real ( kind = rk ) xrite

  f(1:nu) = 0.0D+00
  adiag(1:nu) = 0.0D+00
  aleft(1:nu) = 0.0D+00
  arite(1:nu) = 0.0D+00
!
!  For element IE...
!
  do ie = 1, n

    he = h(ie)
    xleft = xn(node(1,ie))
    xrite = xn(node(2,ie))
!
!  For quadrature point IQ...
!
    do iq = 1, nquad

      xqe = xquad(ie)
!
!  Compute value of U for previous solution.
!
      total = 0.0D+00

      do il = 1, nl

        ig = node(il,ie)
        iu = indx(ig)

        if ( iu <= 0 ) then

          if ( il == 1 ) then
            total = total + ul
          else
            total = total + ur
          end if

        else
          total = total + u(iu)
        end if

      end do

      uold = total / real ( nl, kind = rk )
!
!  Compute value of U' for previous solution.
!
      jl = node(1,ie)
      jr = node(2,ie)
      iul = indx(jl)
      iur = indx(jr)

      if ( iul <= 0 ) then
        uoldx = ( u(iur) - ul ) / he
      else if ( iur <= 0 ) then
        uoldx = ( ur - u(iul) ) / he
      else
        uoldx = ( u(iur) - u(iul) ) / he
      end if
!
!  For basis function IL...
!
      do il = 1, nl

        ig = node(il,ie)
        iu = indx(ig)

        if ( 0 < iu ) then

          call phi ( il, xqe, phii, phiix, xleft, xrite )

          f(iu) = f(iu) + he * phii * ( ff ( xqe, problem ) + uold * uoldx )
!
!  Handle boundary conditions that prescribe the value of U'.
!
          if ( ig == 0 ) then

            x = 0.0D+00
            f(iu) = f(iu) - pp ( problem ) * ul

          else if ( ig == n ) then

            x = 1.0D+00
            f(iu) = f(iu) + pp ( problem ) * ur

          end if
!
!  For basis function JL...
!
          do jl = 1, nl

            jg = node(jl,ie)
            ju = indx(jg)

            call phi ( jl, xqe, phij, phijx, xleft, xrite )

            aij = he * ( pp ( problem ) * phiix * phijx &
                  + qq ( problem ) * phii * phij &
                  + uold * phii * phijx &
                  + uoldx * phij * phii )

            if ( ju <= 0 ) then

              if ( jg == 0 ) then
                f(iu) = f(iu) - aij * ul
              else if ( jg == n ) then
                f(iu) = f(iu) - aij * ur
              end if

            else if ( iu == ju ) then
              adiag(iu) = adiag(iu) + aij
            else if ( ju < iu ) then
              aleft(iu) = aleft(iu) + aij
            else
              arite(iu) = arite(iu) + aij
            end if

          end do

        end if

      end do

    end do

  end do

  return
end
subroutine assemble_picard ( adiag, aleft, arite, f, h, indx, n, nl, &
  node, nquad, nu, problem, u, ul, ur, xn, xquad )

!*****************************************************************************80
!
!! ASSEMBLE_PICARD assembles the Picard linear system.
!
!  Discussion:
!
!    The equation we are trying to solve has the form:
!
!      -d/dx ( p(x) du/dx ) + q(x) * u + u * u' = f(x)
!
!    For the Picard iteration, we need to modify the nonlinear term u * u'
!    so that it is linear in the unknown u, and any other factors of u are
!    lagged.  One way to do this gives us the following equation:
!
!      -d/dx ( p(x) du/dx ) + q(x) * u + u * uold' = f(x)
!
!    where uold is the previous iterate.
!
!    Now we can formulate this system as a (linear) finite element problem
!
!      A * u = rhs
!
!    to be solved for the new approximate solution u.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = rk ) ADIAG(NU).
!    ADIAG(I) is the "diagonal" coefficient of the I-th
!    equation in the linear system.  That is, ADIAG(I) is
!    the coefficient of the I-th unknown in the I-th equation.
!
!    Output, real ( kind = rk ) ALEFT(NU).
!    ALEFT(I) is the "left hand" coefficient of the I-th
!    equation in the linear system.  That is, ALEFT(I) is the
!    coefficient of the (I-1)-th unknown in the I-th equation.
!    There is no value in ALEFT(1), since the first equation
!    does not refer to a "0-th" unknown.
!
!    Output, real ( kind = rk ) ARITE(NU).
!    ARITE(I) is the "right hand" coefficient of the I-th
!    equation in the linear system.  ARITE(I) is the coefficient
!    of the (I+1)-th unknown in the I-th equation.  There is
!    no value in ARITE(NU) because the NU-th equation does not
!    refer to an "NU+1"-th unknown.
!
!    Output, real ( kind = rk ) F(NU), the right hand side of the linear
!    equations.
!
!    Input, real ( kind = rk ) H(N), the length of the subintervals.
!
!    Input, integer INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    Input, integer N.
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    Input, integer NL.
!    The number of basis functions used in a single
!    subinterval.  (NL-1) is the degree of the polynomials
!    used.  For this code, NL is fixed at 2, meaning that
!    piecewise linear functions are used as the basis.
!
!    Input, integer NODE(NL,N).
!    For each subinterval I:
!    NODE(1,I) is the number of the left node, and
!    NODE(2,I) is the number of the right node.
!
!    Input, integer NQUAD.
!    The number of quadrature points used in a subinterval.
!    This code uses NQUAD = 1.
!
!    Input, integer NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Input, real ( kind = rk ) U(NU), the solution value
!    from the previous iteration,
!
!    Input, real ( kind = rk ) UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    Input, real ( kind = rk ) UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    Input, real ( kind = rk ) XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
!    Input, real ( kind = rk ) XQUAD(N)
!    XQUAD(I) is the location of the single quadrature point
!    in interval I.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer nu

  real ( kind = rk ) aij
  real ( kind = rk ) adiag(nu)
  real ( kind = rk ) aleft(nu)
  real ( kind = rk ) arite(nu)
  real ( kind = rk ) f(nu)
  real ( kind = rk ) ff
  real ( kind = rk ) h(n)
  real ( kind = rk ) he
  integer ie
  integer ig
  integer il
  integer indx(0:n)
  integer iu
  integer iul
  integer iur
  integer iq
  integer jg
  integer jl
  integer jr
  integer ju
  integer nl
  integer node(nl,n)
  integer nquad
  real ( kind = rk ) phii
  real ( kind = rk ) phiix
  real ( kind = rk ) phij
  real ( kind = rk ) phijx
  real ( kind = rk ) pp
  integer problem
  real ( kind = rk ) qq
  real ( kind = rk ) total
  real ( kind = rk ) u(nu)
  real ( kind = rk ) ul
  real ( kind = rk ) uold
  real ( kind = rk ) uoldx
  real ( kind = rk ) ur
  real ( kind = rk ) x
  real ( kind = rk ) xleft
  real ( kind = rk ) xn(0:n)
  real ( kind = rk ) xqe
  real ( kind = rk ) xquad(n)
  real ( kind = rk ) xrite

  f(1:nu) = 0.0D+00
  adiag(1:nu) = 0.0D+00
  aleft(1:nu) = 0.0D+00
  arite(1:nu) = 0.0D+00
!
!  For element IE...
!
  do ie = 1, n

    he = h(ie)
    xleft = xn(node(1,ie))
    xrite = xn(node(2,ie))
!
!  For quadrature point IQ...
!
    do iq = 1, nquad

      xqe = xquad(ie)
!
!  Compute value of U for previous solution.
!
      total = 0.0D+00

      do il = 1, nl

        ig = node(il,ie)
        iu = indx(ig)

        if ( iu <= 0 ) then

          if ( il == 1 ) then
            total = total + ul
          else
            total = total + ur
          end if

        else
          total = total + u(iu)
        end if

      end do

      uold = total / real ( nl, kind = rk )
!
!  Compute value of U' for previous solution.
!
      jl = node(1,ie)
      jr = node(2,ie)
      iul = indx(jl)
      iur = indx(jr)

      if ( iul <= 0 ) then
        uoldx = ( u(iur) - ul ) / he
      else if ( iur <= 0 ) then
        uoldx = ( ur - u(iul) ) / he
      else
        uoldx = ( u(iur) - u(iul) ) / he
      end if
!
!  For basis function IL...
!
      do il = 1, nl

        ig = node(il,ie)
        iu = indx(ig)

        if ( 0 < iu ) then

          call phi ( il, xqe, phii, phiix, xleft, xrite )

          f(iu) = f(iu) + he * phii * ( ff ( xqe, problem ) )
!
!  Handle boundary conditions that prescribe the value of U'.
!
          if ( ig == 0 ) then

            x = 0.0
            f(iu) = f(iu) - pp ( problem ) * ul

          else if ( ig == n ) then

            x = 1.0D+00
            f(iu) = f(iu) + pp ( problem ) * ur

          end if
!
!  For basis function JL...
!
          do jl = 1, nl

            jg = node(jl,ie)
            ju = indx(jg)

            call phi ( jl, xqe, phij, phijx, xleft, xrite )

            aij = he * ( pp ( problem ) * phiix * phijx &
                  + qq ( problem ) * phii * phij &
                  + uold * phii * phijx )

            if ( ju <= 0 ) then

              if ( jg == 0 ) then
                f(iu) = f(iu) - aij * ul
              else if ( jg == n ) then
                f(iu) = f(iu) - aij * ur
              end if

            else if ( iu == ju ) then
              adiag(iu) = adiag(iu) + aij
            else if ( ju < iu ) then
              aleft(iu) = aleft(iu) + aij
            else
              arite(iu) = arite(iu) + aij
            end if

          end do

        end if

      end do

    end do

  end do

  return
end
subroutine compare ( u, indx, n, nl, node, nprint, nu, problem, ul, ur, &
  xl, xn, xr )

!*****************************************************************************80
!
!! COMPARE compares the computed and exact solutions.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) U(NU), the solution of the linear equations.
!
!    Input, integer INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    Input, integer N.
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    Input, integer NL.
!    The number of basis functions used in a single
!    subinterval.  (NL-1) is the degree of the polynomials
!    used.  For this code, NL is fixed at 2, meaning that
!    piecewise linear functions are used as the basis.
!
!    Input, integer NODE(NL,N).
!    For each subinterval I:
!    NODE(1,I) is the number of the left node, and
!    NODE(2,I) is the number of the right node.
!
!    Input, integer NPRINT.
!    The number of points at which the computed solution
!    should be printed out when compared to the exact solution.
!
!    Input, integer NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Input, real ( kind = rk ) UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    Input, real ( kind = rk ) UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    Input, real ( kind = rk ) XL.
!    XL is the left endpoint of the interval over which the
!    differential equation is being solved.
!
!    Input, real ( kind = rk ) XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
!    Input, real ( kind = rk ) XR.
!    XR is the right endpoint of the interval over which the
!    differential equation is being solved.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer nl
  integer nu

  real ( kind = rk ) u(nu)
  integer i
  integer ig
  integer indx(0:n)
  integer iu
  integer j
  integer k
  integer node(nl,n)
  integer nprint
  real ( kind = rk ) phii
  real ( kind = rk ) phiix
  integer problem
  real ( kind = rk ) u_exact
  real ( kind = rk ) u_value
  real ( kind = rk ) ul
  real ( kind = rk ) ur
  real ( kind = rk ) ux
  real ( kind = rk ) x
  real ( kind = rk ) xl
  real ( kind = rk ) xleft
  real ( kind = rk ) xn(0:n)
  real ( kind = rk ) xr
  real ( kind = rk ) xrite

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Compare computed and exact solutions:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      X      Computed U      Exact U'
  write ( *, '(a)' ) ' '

  do i = 1, nprint

    x = ( real ( nprint - i,     kind = rk ) * xl   &
        + real (          i - 1, kind = rk ) * xr ) &
        / real ( nprint     - 1, kind = rk )

    ux = u_exact ( x, problem )

    do j = 1, n

      xleft = xn(j-1)
      xrite = xn(j)
!
!  Search for the interval that X lies in.
!
      if ( xleft <= x .and. x <= xrite ) then

        u_value = 0.0D+00

        do k = 1, nl

          ig = node(k,j)
          iu = indx(ig)
          call phi ( k, x, phii, phiix, xleft, xrite )

          if ( iu <= 0 ) then
            if ( j == 1 .and. k == 1 ) then
              u_value = u_value + ul * phii
            else if ( j == n .and. k == nl ) then
              u_value = u_value + ur * phii
            end if
          else
            u_value = u_value + u(iu) * phii
          end if

        end do

        exit

      end if

    end do

    write(*,'(3g14.6)') x, u_value, ux

  end do

  return
end
function ff ( x, problem )

!*****************************************************************************80
!
!! FF returns the right hand side of the differential equation.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, the evaluation point.
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Output, real ( kind = rk ) FF, the value of F(X).
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) ff
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  integer problem
  real ( kind = rk ) x
!
!  Test problem 1
!
  if ( problem == 1 ) then

    ff = x
!
!  Test problem 2
!
  else if ( problem == 2 ) then

    ff = - 0.5D+00 * pi * cos ( 0.5D+00 * pi * x ) &
      + 2.0D+00 * sin ( 0.5D+00 * pi * x ) &
      * ( 1.0D+00 - cos ( 0.5D+00 * pi * x ) ) / pi

  end if

  return
end
subroutine geometry ( h, ibc, indx, nl, node, n, nu, xl, xn, xquad, xr )

!*****************************************************************************80
!
!! GEOMETRY sets up the geometry for the interval [XL,XR].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = rk ) H(N), the length of the subintervals.
!
!    Input, integer IBC.
!    IBC declares what the boundary conditions are.
!    1, at the left endpoint, U has the value UL,
!       at the right endpoint, U' has the value UR.
!    2, at the left endpoint, U' has the value UL,
!       at the right endpoint, U has the value UR.
!    3, at the left endpoint, U has the value UL,
!       and at the right endpoint, U has the value UR.
!    4, at the left endpoint, U' has the value UL,
!       at the right endpoint U' has the value UR.
!
!    Output, integer INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    Input, integer NL.
!    The number of basis functions used in a single
!    subinterval.  (NL-1) is the degree of the polynomials
!    used.  For this code, NL is fixed at 2, meaning that
!    piecewise linear functions are used as the basis.
!
!    Output, integer NODE(NL,N).
!    For each subinterval I:
!    NODE(1,I) is the number of the left node, and
!    NODE(2,I) is the number of the right node.
!
!    Input, integer N.
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    Output, integer NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    Input, real ( kind = rk ) XL.
!    XL is the left endpoint of the interval over which the
!    differential equation is being solved.
!
!    Output, real ( kind = rk ) XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
!    Output, real ( kind = rk ) XQUAD(N)
!    XQUAD(I) is the location of the single quadrature point
!    in interval I.
!
!    Input, real ( kind = rk ) XR.
!    XR is the right endpoint of the interval over which the
!    differential equation is being solved.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nl
  integer n

  real ( kind = rk ) h(n)
  integer i
  integer ibc
  integer indx(0:n)
  integer node(nl,n)
  integer nu
  real ( kind = rk ) xl
  real ( kind = rk ) xn(0:n)
  real ( kind = rk ) xquad(n)
  real ( kind = rk ) xr
!
!  Set the value of XN, the locations of the nodes.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Node      Location'
  write ( *, '(a)' ) ' '
  do i = 0, n
    xn(i) =  ( real ( n - i, kind = rk ) * xl   &
             + real (     i, kind = rk ) * xr ) &
             / real ( n,     kind = rk )
    write(*,'(i6,g14.6)') i, xn(i)
  end do
!
!  Set the lengths of each subinterval.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Subint    Length'
  write ( *, '(a)' ) ' '
  do i = 1, n
    h(i) = xn(i) - xn(i-1)
    write(*,'(i6,g14.6)') i, h(i)
  end do
!
!  Set the quadrature points, each of which is the midpoint of its subinterval.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Subint    Quadrature point'
  write ( *, '(a)' ) ' '
  do i = 1, n
    xquad(i) = 0.5D+00 * ( xn(i-1) + xn(i) )
    write(*,'(i6,g14.6)') i, xquad(i)
  end do
!
!  Set the value of NODE, which records, for each interval,
!  the node numbers at the left and right.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Subint  Left Node  Right Node'
  write ( *, '(a)' ) ' '
  do i = 1, n
    node(1,i) = i - 1
    node(2,i) = i
    write(*,'(3i6)') i, node(1,i), node(2,i)
  end do
!
!  Starting with node 0, see if an unknown is associated with
!  the node.  If so, give it an index.
!
  nu = 0
!
!  Handle first node.
!
  i = 0
  if ( ibc == 1 .or. ibc == 3 ) then
    indx(i) = -1
  else
    nu = nu + 1
    indx(i) = nu
  end if
!
!  Handle nodes 1 through nsub-1
!
  do i = 1, n - 1
    nu = nu + 1
    indx(i) = nu
  end do
!
!  Handle the last node.
!
  i = n
  if ( ibc == 2 .or. ibc == 3 ) then
    indx(i) = -1
  else
    nu = nu + 1
    indx(i) = nu
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Node  Unknown'
  write ( *, '(a)' ) ' '
  do i = 0, n
    write(*,'(2i6)') i, indx(i)
  end do

  return
end
subroutine init ( ibc, imax, nprint, nquad, problem, ul, ur, xl, xr )

!*****************************************************************************80
!
!! INIT initializes variables that define the problem.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Output, integer IBC.
!    IBC declares what the boundary conditions are.
!    1, at the left endpoint, U has the value UL,
!       at the right endpoint, U' has the value UR.
!    2, at the left endpoint, U' has the value UL,
!       at the right endpoint, U has the value UR.
!    3, at the left endpoint, U has the value UL,
!       and at the right endpoint, U has the value UR.
!    4, at the left endpoint, U' has the value UL,
!       at the right endpoint U' has the value UR.
!
!    Output, integer IMAX.
!    The number of Newton iterations to carry out.
!
!    Output, integer NPRINT.
!    The number of points at which the computed solution
!    should be printed out when compared to the exact solution.
!
!    Output, integer NQUAD.
!    The number of quadrature points used in a subinterval.
!    This code uses NQUAD = 1.
!
!    Output, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Output, real ( kind = rk ) UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    Output, real ( kind = rk ) UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    Output, real ( kind = rk ) XL.
!    XL is the left endpoint of the interval over which the
!    differential equation is being solved.
!
!    Output, real ( kind = rk ) XR.
!    XR is the right endpoint of the interval over which the
!    differential equation is being solved.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer ibc
  integer imax
  integer nprint
  integer nquad
  integer problem
  real ( kind = rk ) ul
  real ( kind = rk ) ur
  real ( kind = rk ) xl
  real ( kind = rk ) xr

  ibc = 1
  imax = 10
  nprint = 9
  nquad = 1
  ul = 0.0D+00
  ur = 1.0D+00
  xl = 0.0D+00
  xr = 1.0D+00
!
!  Print out the values that have been set.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'The equation is to be solved for'
  write ( *, '(a,g14.6)' ) 'X greater than XL  =  ', xl
  write ( *, '(a,g14.6)' ) ' and less than XR = ', xr
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'The boundary conditions are:'
  write ( *, '(a)' ) ' '

  if ( ibc == 1 .or. ibc == 3 ) then
    write ( *, '(a,g14.6)' ) '  At X = XL, U=', ul
  else
    write ( *, '(a,g14.6)' ) '  At X = XL, U''=', ul
  end if

  if ( ibc == 2 .or. ibc == 3 ) then
    write ( *, '(a,g14.6)' ) '  At X = XR, U=', ur
  else
    write ( *, '(a,g14.6)' ) '  At X = XR, U''=', ur
  end if

  if ( problem == 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  This is test problem #1:'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  P(X) = 1, Q(X) = 0, F(X) = X.'
    write ( *, '(a)' ) '  Boundary conditions: U(0) = 0, U''(1) = 1.'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  The exact solution is U(X) = X'
  else if ( problem == 2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  This is test problem #2:'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  P(X) = 1, Q(X) = 0, '
    write ( *, '(a)' ) '  F(X) = -0.5*pi*cos(0.5*pi*X)'
    write ( *, '(a)' ) '        + 2*sin(0.5*pi*X)*(1-cos(0.5*pi*X)/pi.'
    write ( *, '(a)' ) '  Boundary conditions: U(0) = 0, U''(1) = 1.'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  The exact solution is U(X) = 2*(1-cos(pi*x/2))/pi'
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) 'Number of quadrature points per element is ', nquad
  write ( *, '(a,i8)' ) 'Number of iterations is ', imax

  return
end
subroutine output ( u, ibc, indx, nsub, nu, ul, ur, xn )

!*****************************************************************************80
!
!! OUTPUT prints out the computed solution at the nodes.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) U(NU), the solution of the linear equations.
!
!    Input, integer IBC.
!    IBC declares what the boundary conditions are.
!    1, at the left endpoint, U has the value UL,
!       at the right endpoint, U' has the value UR.
!    2, at the left endpoint, U' has the value UL,
!       at the right endpoint, U has the value UR.
!    3, at the left endpoint, U has the value UL,
!       and at the right endpoint, U has the value UR.
!    4, at the left endpoint, U' has the value UL,
!       at the right endpoint U' has the value UR.
!
!    Input, integer INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    integer NSUB.
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    Input, integer NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    Input, real ( kind = rk ) UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    Input, real ( kind = rk ) UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    Input, real ( kind = rk ) XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nsub
  integer nu

  integer i
  integer ibc
  integer indx(0:nsub)
  real ( kind = rk ) u(nu)
  real ( kind = rk ) u_value
  real ( kind = rk ) ul
  real ( kind = rk ) ur
  real ( kind = rk ) xn(0:nsub)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Computed solution:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Node    X(I)        U(X(I))'
  write ( *, '(a)' ) ' '

  do i = 0, nsub

    if ( i == 0 ) then

      if ( ibc == 1 .or. ibc == 3 ) then
        u_value = ul
      else
        u_value = u(indx(i))
      end if

    else if ( i == nsub ) then

      if ( ibc == 2 .or. ibc == 3 ) then
        u_value = ur
      else
        u_value = u(indx(i))
      end if

    else

      u_value = u(indx(i))

    end if

    write(*,'(i4,2g14.6)')i, xn(i), u_value

  end do

  return
end
subroutine phi ( il, x, phii, phiix, xleft, xrite )

!*****************************************************************************80
!
!! PHI evaluates a linear basis function and its derivative.
!
!  Discussion:
!
!    In any interval, there are just two basis functions.  The first
!    basis function is a line which is 1 at the left endpoint
!    and 0 at the right.  The second basis function is 0 at
!    the left endpoint and 1 at the right.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, integer IL, the index of the basis function.
!    1, the function which is 1 at XLEFT and 0 at XRITE.
!    2, the function which is 0 at XLEFT and 1 at XRITE.
!
!    Input, real ( kind = rk ) X, the evaluation point.
!
!    Output, real ( kind = rk ) PHII, PHIIX, the value of the
!    basis function and its derivative at X.
!
!    Input, real ( kind = rk ) XLEFT, XRITE, the left and right
!    endpoints of the interval.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer il
  real ( kind = rk ) phii
  real ( kind = rk ) phiix
  real ( kind = rk ) x
  real ( kind = rk ) xleft
  real ( kind = rk ) xrite

  if ( xleft <= x .and. x <= xrite ) then

    if ( il == 1 ) then
      phii = ( xrite - x ) / ( xrite - xleft )
      phiix = -1.0D+00 / ( xrite - xleft )
    else
      phii = ( x - xleft ) / ( xrite - xleft )
      phiix = 1.0D+00 / ( xrite - xleft )
    end if
!
!  If X is outside of the interval, then the basis function
!  is always zero.
!
  else

    phii = 0.0D+00
    phiix = 0.0D+00

  end if

  return
end
function pp ( problem )

!*****************************************************************************80
!
!! PP returns the value of the coefficient function P(X).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, 
!      u = x, 
!      p=1, 
!      q=0,
!      f=x,
!      u(0)=0,
!      u'(1)=1.
!
!    * 2, 
!      u = 2*(1-cos(0.5*pi*x))/pi, 
!      p=1, 
!      q=0,
!      f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!      u(0) = 0,
!      u'(1)=1.
!
!    Output, real ( kind = rk ) PP, the value of P(X).
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) pp
  integer problem
!
!  Test problem 1
!
  if ( problem == 1 ) then

    pp = 1.0D+00
!
!  Test problem 2
!
  else if ( problem == 2 ) then

    pp = 1.0D+00

  end if

  return
end
subroutine prsys ( adiag, aleft, arite, f, nu )

!*****************************************************************************80
!
!! PRSYS prints out the tridiagonal linear system to be solved.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) ADIAG(NU), ALEFT(NU), ARITE(NU),
!    the diagonal, left and right entries of the equations.
!
!    Input, real ( kind = rk ) F(NU), the right hand side of the linear
!    system to be solved.
!
!    Input, integer NU.
!    NU is the number of equations to be solved.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nu

  real ( kind = rk ) adiag(nu)
  real ( kind = rk ) aleft(nu)
  real ( kind = rk ) arite(nu)
  real ( kind = rk ) f(nu)
  integer i

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Printout of tridiagonal linear system:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Equation  ALEFT  ADIAG  ARITE  RHS'
  write ( *, '(a)' ) ' '
  do i = 1, nu
    if ( i == 1 ) then
      write(*,'(i3,14x,3g14.6)') i, adiag(i), arite(i), f(i)
    else if ( i < nu ) then
      write(*,'(i3,4g14.6)') i, aleft(i), adiag(i), arite(i), f(i)
    else
      write(*,'(i3,2g14.6,14x,g14.6)') i, aleft(i), adiag(i), f(i)
    end if
  end do

  return
end
function qq ( problem )

!*****************************************************************************80
!
!! QQ returns the value of the coefficient function Q(X).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Output, real ( kind = rk ) QQ, the value of Q(X).
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer problem
  real ( kind = rk ) qq
!
!  Test problem 1
!
  if ( problem == 1 ) then

    qq = 0.0D+00
!
!  Test problem 2
!
  else if ( problem == 2 ) then

    qq = 0.0D+00

  end if

  return
end
subroutine solve ( adiag, aleft, arite, f, nu, u )

!*****************************************************************************80
!
!! SOLVE solves a tridiagonal matrix system of the form A*x = b.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input/output, real ( kind = rk ) ADIAG(NU), ALEFT(NU), ARITE(NU).
!    On input, ADIAG, ALEFT, and ARITE contain the diagonal,
!    left and right entries of the equations.
!    On output, ADIAG and ARITE have been changed in order
!    to compute the solution.
!    Note that for the first equation, there is no ALEFT
!    coefficient, and for the last, there is no ARITE.
!    So there is no need to store a value in ALEFT(1), nor
!    in ARITE(NU).
!
!    Input, real ( kind = rk ) F(NU), the right hand side of the linear
!    system to be solved.
!
!    Input, integer NU.
!    NU is the number of equations to be solved.
!
!    Output, real ( kind = rk ) U(NU), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer nu

  real ( kind = rk ) adiag(nu)
  real ( kind = rk ) aleft(nu)
  real ( kind = rk ) arite(nu-1)
  real ( kind = rk ) f(nu)
  integer i
  real ( kind = rk ) u(nu)

  arite(1) = arite(1) / adiag(1)
  do i = 2, nu-1
    adiag(i) = adiag(i) - aleft(i) * arite(i-1)
    arite(i) = arite(i) / adiag(i)
  end do
  adiag(nu) = adiag(nu) - aleft(nu) * arite(nu-1)

  u(1:nu) = f(1:nu)

  u(1) = u(1) / adiag(1)
  do i = 2, nu
    u(i) = ( u(i) - aleft(i) * u(i-1) ) / adiag(i)
  end do

  do i = nu-1, 1, -1
    u(i) = u(i) - arite(i) * u(i+1)
  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
function u_exact ( x, problem )

!*****************************************************************************80
!
!! U_EXACT returns the value of the exact solution at a point X.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 November 2006
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, the evaluation point.
!
!    Input, integer PROBLEM, indicates which problem to be solved.
!    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
!    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
!    u(0) = 0, u'(1)=1.
!
!    Output, real ( kind = rk ) U_EXACT, the value of the exact solution at X.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  integer problem
  real ( kind = rk ) u_exact
  real ( kind = rk ) x
!
!  Test problem 1
!
  if ( problem == 1 ) then

    u_exact = x
!
!  Test problem 2
!
  else if ( problem == 2 ) then

    u_exact = 2.0D+00 * ( 1.0D+00 - cos ( 0.5D+00 * pi * x ) ) / pi

  end if

  return
end