program onedfe
!
!*******************************************************************************
!
!! ONEDFE solves a one dimensional ODE using finite element methods.
!
!  Discussion:
!
!    The differential equation solved is
!
!      - d/dX (P dU/dX) + Q U  =  F
!
!    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.
!
!
!    The interval [XL,XR] is "meshed" with NSUB+1 points,
!
!    XN(0) = XL, XN(1)=XL+H, XN(2)=XL+2*H, ..., XN(NSUB)=XR.
!
!    This creates NSUB subintervals, with interval number 1
!    having endpoints XN(0) and XN(1), and so on up to interval
!    NSUB, which has endpoints XN(NSUB-1) and XN(NSUB).
!
!  Parameters:
!
!    REAL ADIAG(NU), the "diagonal" coefficients.  That is, ADIAG(I) is
!    the coefficient of the I-th unknown in the I-th equation.
!
!    REAL ALEFT(NU), the "left hand" coefficients.  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 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 F(NSUB+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.
!
!  REAL H(NSUB)
!     H(I) is the length of subinterval I.  This code uses
!     equal spacing for all 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 INDX(0:NSUB).
!     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 NSUB, which will be unknown NSUB or NSUB+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(NSUB,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.
!
!  INTEGER NQUAD.
!     The number of quadrature points used in a subinterval.
!     This code uses NQUAD = 1.
!
!  INTEGER NSUB.
!     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 NSUB-1,
!     NSUB, or NSUB+1 unknown values, which are the coefficients
!     of basis functions.
!
!  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.
!
!  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.
!
!  REAL XL.
!     XL is the left endpoint of the interval over which the
!     differential equation is being solved.
!
!  REAL XN(0:NSUB).
!     XN(I) is the location of the I-th node.  XN(0) is XL,
!     and XN(NSUB) is XR.
!
!  REAL XQUAD(NSUB)
!     XQUAD(I) is the location of the single quadrature point
!     in interval I.
!
!  REAL XR.
!     XR is the right endpoint of the interval over which the
!     differential equation is being solved.
!
  implicit none
!
!  Set the number of subintervals.
!
  integer, parameter :: nsub = 5
!
!  Set the number of basis functions per interval.
!  2 means linear functions are used.
!
  integer, parameter :: nl = 2
!
  real adiag(nsub+1)
  real aleft(nsub+1)
  real arite(nsub+1)
  real f(nsub+1)
  real h(nsub)
  integer ibc
  integer indx(0:nsub)
  integer node(nsub,nl)
  integer nquad
  integer nu
  real ul
  real ur
  real xl
  real xn(0:nsub)
  real xquad(nsub)
  real xr
!
  call timestamp ( )

  print *, ' '
  print *, 'ONEDFE'
  print *, ' '
  print *, 'Solve the two-point boundary value problem'
  print *, ' '
  print *, '  - d/dX (P dU/dX) + Q U  =  F'
  print *, ' '
  print *, 'on the interval [XL,XR], specifying'
  print *, 'the value of U or U'' at each end.'
  print *, ' '
  print *, 'The interval [XL,XR] is broken into ', nsub, ' subintervals'
  print *, 'Number of basis functions per element is ',nl
!
!  Initialize the data that defines the problem.
!
  call init ( ibc, nquad, ul, ur, xl, xr )
!
!  Compute the quantities which define the geometry of the
!  problem.
!
  call geometry ( h, ibc, indx, nl, node, nsub, nu, xl, xn, xquad, xr )
!
!  Assemble the linear system.
!
  call assemble ( adiag, aleft, arite, f, h, indx, nl, node, nu, nquad, &
    nsub, ul, ur, xn, xquad )
!
!  Print out the linear system.
!
  call prsys ( adiag, aleft, arite, f, nu )
!
!  Solve the linear system.
!
  call solve ( adiag, aleft, arite, f, nu )
!
!  Print out the solution.
!
  call output ( f, ibc, indx, nsub, nu, ul, ur, xn )

  stop
end
subroutine assemble ( adiag, aleft, arite, f, h, indx, nl, node, nu, &
  nquad, nsub, ul, ur, xn, xquad )
!
!*******************************************************************************
!
!! ASSEMBLE assembles the matrix and right-hand-side of the linear system.
!
!
!  Discussion:
!
!    The linear system has the form:
!
!      K * C  =  F
!
!    that is to be solved for the coefficients C.
!
!    Numerical integration is used to compute the entries of K and F.
!
  implicit none
!
  integer nl
  integer nsub
  integer nu
!
  real aij
  real adiag(nu)
  real aleft(nu)
  real arite(nu)
  real f(nu)
  real ff
  real h(nsub)
  real he
  integer i
  integer ie
  integer ig
  integer il
  integer indx(0:nsub)
  integer iq
  integer iu
  integer jg
  integer jl
  integer ju
  integer node(nsub,nl)
  integer nquad
  real phii
  real phiix
  real phij
  real phijx
  real pp
  real qq
  real ul
  real ur
  real xleft
  real xn(0:nsub)
  real xquad(nsub)
  real xquade
  real xrite
!
!  Zero out the arrays that hold the coefficients of the matrix
!  and the right hand side.
!
  f(1:nu) = 0.0E+00
  adiag(1:nu) = 0.0E+00
  aleft(1:nu) = 0.0
  arite(1:nu) = 0.0
!
!  For interval number IE,
!
  do ie = 1, nsub

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

      xquade = xquad(ie)
!
!  and evaluate the integrals associated with the basis functions
!  for the left, and for the right nodes.
!
      do il = 1, nl

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

        if ( iu > 0 ) then

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

          f(iu) = f(iu) + he * ff(xquade) * phii
!
!  Take care of boundary nodes at which U' was specified.
!
          if ( ig == 0 ) then

            f(iu) = f(iu) - pp(0.0E+00)*ul

          else if ( ig == nsub ) then

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

          end if
!
!  Evaluate the integrals that take a product of the basis
!  function times itself, or times the other basis function
!  that is nonzero in this interval.
!
          do jl = 1, nl

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

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

            aij = he * ( pp(xquade)*phiix*phijx + qq(xquade)*phii*phij   )
!
!  If there is no variable associated with the node, then it's
!  a specified boundary value, so we multiply the coefficient
!  times the specified boundary value and subtract it from the
!  right hand side.
!
            if ( ju <= 0 ) then

              if ( jg == 0 ) then

                f(iu) = f(iu) - aij * ul

              else if ( jg == nsub ) then

                f(iu) = f(iu) - aij * ur

              end if
!
!  Otherwise, we add the coefficient we've just computed to the
!  diagonal, or left or right entries of row IU of the matrix.
!
            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 if

          end do

        end if

      end do
    end do
  end do

  return
end
function ff ( x )
!
!*******************************************************************************
!
!! FF evaluates the right hand side function.
!
!
!  Discussion:
!
!    This routine evaluates the function F(X) in the differential equation.
!
!      -d/dx (p du/dx) + q u  =  f
!
!    at the point X.
!
!  Parameters:
!
!    Input, real X, the argument of the function.
!
!    Output, real FF, the value of the function.
!
  implicit none
!
  real ff
  real x
!
  ff = 0.0E+00

  return
end
subroutine geometry ( h, ibc, indx, nl, node, nsub, nu, xl, xn, xquad, &
  xr )
!
!*******************************************************************************
!
!! GEOMETRY 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)
    print *, 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)
    print *, 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.5E+00 * (xn(i-1)+xn(i))
    print *, 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
    print *, 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
    print *,  i, indx(i)
  end do

  return
end
subroutine init ( ibc, nquad, ul, ur, xl, xr )
!
!*******************************************************************************
!
!! INIT assigns values to variables which define the problem.
!
  implicit none
!
  integer ibc
  integer nquad
  real ul
  real ur
  real xl
  real xr
!
!  IBC declares what the boundary conditions are.
!
  ibc = 1
!
!  NQUAD is the number of quadrature points per subinterval.
!  The program as currently written cannot handle any value for
!  NQUAD except 1!
!
  nquad = 1
!
!  Set the values of U or U' at the endpoints.
!
  ul = 0.0E+00
  ur = 1.0E+00
!
!  Define the location of the endpoints of the interval.
!
  xl = 0.0E+00
  xr = 1.0E+00
!
!  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

  return
end
subroutine output ( f, ibc, indx, nsub, nu, ul, ur, xn )
!
!*******************************************************************************
!
!! OUTPUT prints out the computed solution.
!
!
!  Discussion:
!
!    We simply print out the solution vector F, except that, for
!    certain boundary conditions, we are going to have to get the
!    value of the solution at XL or XR by using the specified
!    boundary value.
!
  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)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Computed solution coefficients:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Node    X(I)        U(X(I))'
  write ( *, '(a)' ) ' '

  do i = 0, nsub
!
!  If we're at the first node, check the boundary condition.
!
    if ( i == 0 ) then

      if ( ibc == 1 .or. ibc == 3 ) then
        u = ul
      else
        u = f(indx(i))
      end if
!
!  If we're at the last node, check the boundary condition.
!
    else if ( i == nsub ) then

      if ( ibc == 2 .or. ibc == 3 ) then
        u = ur
      else
        u = f(indx(i))
      end if
!
!  Any other node, we're sure the value is stored in F.
!
    else

      u = f(indx(i))

    end if

    write ( *, '(i6,f6.2,g14.6)' ) i, xn(i), u

  end do

  return
end
subroutine phi ( il, x, phii, phiix, xleft, xrite )
!
!*******************************************************************************
!
!! PHI evaluates a linear basis function and its derivative.
!
!
!  The evaluation is done at a point X in an interval [XLEFT,XRITE].
!
!  In this interval, there are just two nonzero 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.0E+00 / (xrite-xleft)
    else
      phii = (x-xleft) / (xrite-xleft)
      phiix = 1.0E+00 / (xrite-xleft)
    end if
!
!  If X is outside of the interval, just set everything to 0.
!
  else

    phii = 0.0E+00
    phiix = 0.0E+00

  end if

  return
end
function pp ( x )
!
!*******************************************************************************
!
!! PP evaluates the function P in the differential equation.
!
!
!  Discussion:
!
!    The function P appears in the differential equation as;
!
!      - d/dx (p du/dx) + q u  =  f
!
!  Parameters:
!
!    Input, real X, the argument of the function.
!
!    Output, real PP, the value of the function.
!
  implicit none
!
  real pp
  real x
!
  pp = 1.0E+00

  return
end
subroutine prsys ( adiag, aleft, arite, f, nu )
!
!*******************************************************************************
!
!! PRSYS prints out the tridiagonal linear system.
!
  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
    print *, i, aleft(i), adiag(i), arite(i), f(i)
  end do

  return
end
function qq ( x )
!
!*******************************************************************************
!
!! QQ evaluates the function Q in the differential equation.
!
!
!  Discussion:
!
!    The function Q appears in the differential equation as:
!
!      - d/dx (p du/dx) + q u  =  f
!
!  Parameters:
!
!    Input, real X, the argument of the function.
!
!    Output, real QQ, the value of the function.
!
  implicit none
!
  real qq
  real x
!
  qq = 0.0E+00

  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, 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
!
!  Carry out Gauss elimination on the matrix, saving information
!  needed for the backsolve.
!
  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)
!
!  Carry out the same elimination steps on F that were done to the
!  matrix.
!
  f(1) = f(1) / adiag(1)
  do i = 2, nu
    f(i) = ( f(i) - aleft(i) * f(i-1) ) / adiag(i)
  end do
!
!  And now carry out the steps of "back substitution".
!
  do i = nu-1, 1, -1
    f(i) = f(i)-arite(i)*f(i+1)
  end do

  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