program nonlin
!
!*******************************************************************************
!
!! NONLIN solves a nonlinear one dimensional boundary value problem.
!
!
!  Discussion:
!
!    The differential equation has the form:
!
!      -d/dx (p(x) du/dx) + q(x)*u  =  f(x) - u*u'
!
!    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 problems:
!
!    u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
!    The code should solve this problem exactly.
!
!    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.
!
!  Parameters:
!
!  ADIAG  REAL 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.
!
!  ALEFT  REAL 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.
!
!  ARITE  REAL 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.
!
!  F  REAL F(N+1) or F(NU).
!
!     ASSEMB stores into F the right hand side of the linear
!     equations.
!
!     SOLVE replaces those values of F by the solution of the
!     linear equations.
!
!  FOLD   REAL FOLD(N+1) or FOLD(NU).
!     FOLD contains the value of F from the previous iteration,
!     and is used in ASSEMB to add correction terms to the
!     matrix and right hand side.
!
!  H  REAL H(N)
!     H(I) is the length of subinterval I.  This code uses
!     equal spacing for all the subintervals.
!
!  IBC    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.
!
!  IMAX   INTEGER IMAX.
!     The number of iterations to carry out.
!
!  INDX   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.
!
!  NL     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.
!
!  NODE   INTEGER NODE(N,NL).
!     For each subinterval I:
!
!       NODE(I,1) is the number of the left node, and
!       NODE(I,2) is the number of the right node.
!
!  NPRINT INTEGER NPRINT.
!     The number of points at which the computed solution
!     should be printed out when compared to the exact
!     solution.
!
!  NQUAD  INTEGER NQUAD.
!     The number of quadrature points used in a subinterval.
!     This code uses NQUAD = 1.
!
!  N  INTEGER N.
!     The number of subintervals into which the interval
!     [XL,XR] is broken.
!
!  NU     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.
!
!  UL     REAL 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.
!
!  UR     REAL 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.
!
!  XL     REAL XL.
!     XL is the left endpoint of the interval over which the
!     differential equation is being solved.
!
!  XN     REAL XN(0:N).
!     XN(I) is the location of the I-th node.  XN(0) is XL,
!     and XN(N) is XR.
!
!  XNEWT  REAL XNEWT.
!     A factor which is 0 for the first three iterations, and
!     1 thereafter.  It multiplies some terms in ASSEMB,
!     and thus is used to "weaken" the nonlinearity for
!     the first few steps.
!
!  XQUAD  REAL XQUAD(N)
!     XQUAD(I) is the location of the single quadrature point
!     in interval I.
!
!  XR     REAL XR.
!     XR is the right endpoint of the interval over which the
!     differential equation is being solved.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: nl = 2
!
  real adiag(n+1)
  real aleft(n+1)
  real arite(n+1)
  real f(n+1)
  real fold(n+1)
  real h(n)
  integer i
  integer ibc
  integer imax
  integer indx(0:n)
  integer j
  integer node(n,nl)
  integer nprint
  integer nquad
  integer nu
  real ul
  real ur
  real xl
  real xn(0:n)
  real xnewt
  real xquad(n)
  real xr
!
  call timestamp ( )

  print *, ' '
  print *, 'NONLIN'
  print *, '  Solve a nonlinear boundary value problem:'
  print *, '    -d/dx (p(x) du/dx) + q(x)*u  =  f(x) - u*u'' '
  print *, '  on an interval [xl,xr], with the values of'
  print *, '  u or u'' specified at xl and xr.'
!
!  Initialize variables that define the problem.
!
  call init ( ibc, imax, nprint, nquad, ul, ur, xl, xr )
!
!  Compute the quantities that describe the geometry of the problem.
!
  call geomet ( h, ibc, indx, nl, node, n, nu, xl, xn, xquad, xr )
!
!  Initialize the "previous" solution to 0.
!
  fold(1:nu) = 0.0E+00
!
!  Begin the iteration.
!
  do i = 1, imax
!
!  Is it time for full nonlinear Newton iteration?
!
    if ( i <= 3 ) then
      xnewt = 0.0
    else
      xnewt = 1.0
    end if
!
!  Assemble the matrix.
!
    call assemb ( adiag, aleft, arite, f, fold, h, indx, n, nl, node, &
      nquad, nu, ul, ur, xn, xnewt, xquad )
!
!  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 )
!
!  Print the current solution.
!
    call output ( f, ibc, indx, n, nu, ul, ur, xn )
!
!  Save a copy of the current solution in FOLD.
!
    fold(1:nu) = f(1:nu)

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

  stop
end
subroutine assemb ( adiag, aleft, arite, f, fold, h, indx, n, nl, &
  node, nquad, nu, ul, ur, xn, xnewt, xquad )
!
!*******************************************************************************
!
!! ASSEMB assembles the matrix and right hand side of the linear system.
!
  implicit none
!
  integer n
  integer nu
!
  real aij
  real adiag(nu)
  real aleft(nu)
  real arite(nu)
  real f(nu)
  real ff
  real fold(nu)
  real h(n)
  real he
  integer i
  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(n,nl)
  integer nquad
  real phii
  real phiix
  real phij
  real phijx
  real pp
  real qq
  real sum
  real ul
  real uold
  real uoldx
  real ur
  real xleft
  real xn(0:n)
  real xnewt
  real xqe
  real xquad(n)
  real xrite
!
  f(1:nu) = 0.0E+00
  adiag(1:nu) = 0.0E+00
  aleft(1:nu) = 0.0E+00
  arite(1:nu) = 0.0E+00
!
!  For element IE...
!
  do ie = 1, n

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

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

      do il = 1, nl

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

        if ( iu <= 0 ) then

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

        else
          sum = sum+fold(iu)
        end if

      end do

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

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

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

        if ( iu > 0 ) then

          call phi ( il, xqe, phii, phiix, xleft, xrite )
          f(iu) = f(iu) + he*phii*(ff(xqe)+uold*uoldx*xnewt)
!
!  Handle boundary conditions that prescribe the value of U'.
!
          if ( ig == 0 ) then

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

          else if ( ig == n ) then

            f(iu) = f(iu) + pp(1.0)*ur

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

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

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

            aij = he*( pp(xqe)*phiix*phijx &
                  +qq(xqe)*phii*phij &
                  +uold*phii*phijx &
                  +xnewt*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 ( iu > ju ) 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
function ff ( x )
!
!*******************************************************************************
!
!! FF returns the right hand side of the differential equation.
!
  implicit none
!
  real ff
  real, parameter :: pi = 3.14159265E+00
  real x
!
!  Test problem 1
!
!     ff = x
!
!  Test problem 2
!
  ff = - 0.5 * pi * cos ( 0.5 * pi * x ) &
    + 2 * sin ( 0.5 * pi * x ) * ( 1.0 - cos ( 0.5 * pi * x ) ) / pi

  return
end
subroutine geomet ( h, ibc, indx, nl, node, nsub, nu, xl, xn, xquad, &
  xr )
!
!*******************************************************************************
!
!! GEOMET sets up the geometry for the interval [XL,XR].
!
  implicit none
!
  integer nl
  integer nsub
!
  real h(nsub)
  integer i
  integer ibc
  integer indx(0:nsub)
  integer node(nsub,nl)
  integer nu
  real xl
  real xn(0:nsub)
  real xquad(nsub)
  real xr
!
!  Set the value of XN, the locations of the nodes.
!
  print *, ' '
  print *, '  Node      Location'
  print *, ' '
  do i = 0, nsub
    xn(i) =  ( (nsub-i)*xl + i*xr ) / real(nsub)
    write(*,'(i6,g14.6)')i,xn(i)
  end do
!
!  Set the lengths of each subinterval.
!
  print *, ' '
  print *, 'Subint    Length'
  print *, ' '
  do i = 1, nsub
    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.
!
  print *, ' '
  print *, 'Subint    Quadrature point'
  print *, ' '
  do i = 1, nsub
    xquad(i) = 0.5*(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.
!
  print *, ' '
  print *, 'Subint  Left Node  Right Node'
  print *, ' '
  do i = 1, nsub
    node(i,1) = i-1
    node(i,2) = i
    write(*,'(3i6)')i,node(i,1),node(i,2)
  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,nsub-1
    nu = nu+1
    indx(i) = nu
  end do
!
!  Handle the last node.
!
  i = nsub
  if ( ibc == 2 .or. ibc == 3 ) then
    indx(i) = -1
   else
    nu = nu+1
    indx(i) = nu
   end if

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

  return
end
subroutine init ( ibc, imax, nprint, nquad, ul, ur, xl, xr )
!
!*******************************************************************************
!
!! INIT initializes variables that define the problem.
!
  implicit none
!
  integer ibc
  integer imax
  integer nprint
  integer nquad
  real ul
  real ur
  real xl
  real xr
!
  ibc = 1
  imax = 5
  nprint = 9
  nquad = 1
  ul = 0.0
  ur = 1.0
  xl = 0.0
  xr = 1.0
!
!  Print out the values that have been set.
!
  print *, ' '
  print *, 'The equation is to be solved for'
  print *, 'X greater than XL  =  ',xl,' and less than XR = ',xr
  print *, ' '
  print *, 'The boundary conditions are:'
  print *, ' '
  if ( ibc == 1 .or. ibc == 3 ) then
    print *, '  At X = XL, U=',ul
  else
    print *, '  At X = XL, U''=',ul
  end if
  if ( ibc == 2 .or. ibc == 3 ) then
    print *, '  At X = XR, U=',ur
  else
    print *, '  At X = XR, U''=',ur
  end if
  print *, ' '
  print *, 'Number of quadrature points per element is ',nquad
  print *, 'Number of iterations is ',imax

  return
end
subroutine output ( f, ibc, indx, nsub, nu, ul, ur, xn )
!
!*******************************************************************************
!
!! OUTPUT prints out the computed solution at the nodes.
!
  implicit none
!
  integer nsub
  integer nu
!
  real f(nu)
  integer i
  integer ibc
  integer indx(0:nsub)
  real u
  real ul
  real ur
  real xn(0:nsub)
!
  print *, ' '
  print *, 'Computed solution:'
  print *, ' '
  print *, 'Node    X(I)        U(X(I))'
  print *, ' '

  do i = 0, nsub
    if ( i == 0 ) then
      if ( ibc == 1 .or. ibc == 3 ) then
        u = ul
      else
        u = f(indx(i))
      end if
    else if ( i == nsub ) then
      if ( ibc == 2 .or. ibc == 3 ) then
        u = ur
      else
        u = f(indx(i))
      end if
    else
      u = f(indx(i))
    end if
    write(*,'(i4,2g14.6)')i,xn(i),u

  end do

  return
end
subroutine out2 ( f, indx, n, nl, node, nprint, nu, ul, ur, xl, xn, &
  xr )
!
!*******************************************************************************
!
!! OUT2 compares the computed and exact solutions.
!
  implicit none
!
  integer n
  integer nl
  integer nu
!
  real f(nu)
  integer i
  integer ig
  integer indx(0:n)
  integer iu
  integer j
  integer k
  integer node(n,nl)
  integer nprint
  real phii
  real phiix
  real u
  real uex
  real ul
  real ur
  real ux
  real x
  real xl
  real xleft
  real xn(0:n)
  real xr
  real xrite
!
  print *, ' '
  print *, 'Compare computed and exact solutions:'
  print *, ' '
  print *, '      X      Computed U      Exact U'
  print *, ' '
  do i = 1,nprint
    x = ( (nprint-i)*xl + (i-1)*xr) / real(nprint-1)
    ux = uex(x)
    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 = 0.0

        do k = 1,nl
          ig = node(j,k)
          iu = indx(ig)
          call phi(k,x,phii,phiix,xleft,xrite)
          if ( iu <= 0 ) then
            if ( j == 1 .and. k == 1 ) then
              u = u+ul*phii
            else if ( j == n .and. k == nl ) then
              u = u+ur*phii
            end if
          else
            u = u+f(iu)*phii
          end if
        end do

        exit

      end if

    end do

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

  end do

  return
end
subroutine phi ( il, x, phii, phiix, xleft, xrite )
!
!*******************************************************************************
!
!! PHI evaluates a linear basis function and its derivative.
!
!
!  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.
!
  implicit none
!
  integer il
  real phii
  real phiix
  real x
  real xleft
  real xrite
!
  if ( xleft <= x .and. x <= xrite ) then
    if ( il == 1 ) then
      phii = (xrite-x)/(xrite-xleft)
      phiix = -1.0/(xrite-xleft)
    else
      phii = (x-xleft)/(xrite-xleft)
      phiix = 1.0/(xrite-xleft)
    end if
!
!  If X is outside of the interval, then the basis function
!  is always zero.
!
  else
    phii = 0.0
    phiix = 0.0
  end if

  return
end
function pp ( x )
!
!*******************************************************************************
!
!! PP returns the value of the coefficient function P(X).
!
  implicit none
!
  real pp
  real x
!
!  Test problem 1
!
  pp = 1.0
!
!  Test problem 2
!
  pp = 1.0

  return
end
subroutine prsys ( adiag, aleft, arite, f, nu )
!
!*******************************************************************************
!
!! PRSYS prints out the tridiagonal linear system to be solved.
!
  implicit none
!
  integer nu
!
  real adiag(nu)
  real aleft(nu)
  real arite(nu)
  real f(nu)
  integer i
!
  print *, ' '
  print *, 'Printout of tridiagonal linear system:'
  print *, ' '
  print *, 'Equation  ALEFT  ADIAG  ARITE  RHS'
  print *, ' '
  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 ( x )
!
!*******************************************************************************
!
!! QQ returns the value of the coefficient function Q(X).
!
  implicit none
!
  real qq
  real x
!
!  Test problem 1
!
  qq = 0.0
!
!  Test problem 2
!
  qq = 0.0

  return
end
subroutine solve ( adiag, aleft, arite, f, nu )
!
!*******************************************************************************
!
!! SOLVE solves a tridiagonal matrix system of the form A*x = b.
!
!
!  Parameters:
!
!    Input/output, REAL 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/output, REAL F(NU).
!    On input, F contains the right hand side of the linear
!    system to be solved.
!    On output, F contains the solution of the linear system.
!
!    Input, INTEGER NU.
!    NU is the number of equations to be solved.
!
  implicit none
!
  integer nu
!
  real adiag(nu)
  real aleft(nu)
  real arite(nu-1)
  real f(nu)
  integer i
!
  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)

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

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

  return
end
function uex ( x )
!
!*******************************************************************************
!
!! UEX returns the value of the exact solution at a point X.
!
  implicit none
!
  real, parameter :: pi = 3.14159265E+00
  real uex
  real x
!
!  Test problem 1
!
!     uex = x
!
!  Test problem 2
!
  uex = 2*(1-cos(0.5*pi*x))/pi

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  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
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end