difequ Subroutine

subroutine difequ ( mode, xx, v )

  !*****************************************************************************80
  !
  !! DIFEQU returns information about a differential equation.
  !
  !  Discussion:
  !
  !    This sample version of DIFEQU is for the example in chapter XV.  It is a
  !    nonlinear second order two point boundary value problem.
  !
  !  Modified:
  !
  !    16 February 2007
  !
  !  Author:
  !
  !    Carl DeBoor
  !
  !  Reference:
  !
  !    Carl DeBoor,
  !    A Practical Guide to Splines,
  !    Springer, 2001,
  !    ISBN: 0387953663,
  !    LC: QA1.A647.v27.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) MODE, an integer indicating the task to 
  !    be performed.
  !    1, initialization
  !    2, evaluate the differential equation at point XX.
  !    3, specify the next side condition
  !    4, analyze the approximation
  !
  !    Input, real ( kind = 8 ) XX, a point at which information is wanted
  !
  !    Output, real ( kind = 8 ) V, depends on the MODE.
  !
  implicit none

  integer ( kind = 4 ), parameter :: npiece = 100
  integer ( kind = 4 ), parameter :: ncoef = 2000

  real ( kind = 8 ) break
  real ( kind = 8 ) coef
  real ( kind = 8 ), save :: eps
  real ( kind = 8 ) ep1
  real ( kind = 8 ) ep2
  real ( kind = 8 ) error
  real ( kind = 8 ), save :: factor
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iside
  integer ( kind = 4 ) itermx
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kpm
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mode
  real ( kind = 8 ) ppvalu
  real ( kind = 8 ) rho
  real ( kind = 8 ), save :: s2ovep
  real ( kind = 8 ) solutn
  real ( kind = 8 ) un
  real ( kind = 8 ) v(20)
  real ( kind = 8 ) value
  real ( kind = 8 ) x
  real ( kind = 8 ) xside
  real ( kind = 8 ) xx

  save / approx /
  save / other /
  save / side /

  common / approx / break(npiece), coef(ncoef), l, kpm
  common / other / itermx, k, rho(19)
  common / side / m, iside, xside(10)
  !
  !  Initialize everything,  Set the order M of the differential equation, 
  !  the nondecreasing sequence XSIDE(1:M), of points at which side 
  !  conditions are given and anything else necessary.
  !
  if ( mode == 1 ) then

     m = 2
     xside(1) = 0.0D+00
     xside(2) = 1.0D+00
     !
     !  Print out heading.
     !
     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) '  Carrier''s nonlinear perturbation problem'
     write ( *, '(a)' ) ' '

     eps = 0.005D+00
     write ( *, '(a,g14.6)' ) '  EPS = ', eps
     !
     !  Set constants used in formula for solution below.
     !
     factor = ( sqrt ( 2.0D+00 ) + sqrt ( 3.0D+00 ) )**2
     s2ovep = sqrt ( 2.0D+00 / eps )
     !
     !  Initial guess for Newton iteration: UN(X) = X*X-1.
     !
     l = 1
     break(1) = 0.0D+00
     coef(1:kpm) = 0.0D+00
     coef(1) = -1.0D+00
     coef(3) = 2.0D+00
     itermx = 10
     !
     !  Provide value of left side coefficients and right hand side at XX.
     !  Specifically, at XX the differential equation reads:
     !
     !    V(M+1) D**M + V(M) D**(M-1) + ... + V(1) D**0 = V(M+2)
     !
     !  in terms of the quantities V(1:M+2), to be computed here.
     !
  else if ( mode == 2 ) then

     v(3) = eps
     v(2) = 0.0D+00

     un = ppvalu ( break, coef, l, kpm, xx, 0 )

     v(1) = 2.0D+00 * un
     v(4) = un**2 + 1.0D+00
     !
     !  Provide the M side conditions. these conditions are of the form
     !
     !    V(M+1) D**M + V(M) D**(M-1) + ... + V(1) D**0 = V(M+2)
     !
     !  in terms of the quantities V(1:M+2), to be specified here.
     !  Note that V(M+1) = 0 for customary side conditions.
     !
  else if ( mode == 3 ) then

     v(m+1) = 0.0D+00

     if ( iside == 1 ) then
        v(2) = 1.0D+00
        v(1) = 0.0D+00
        v(4) = 0.0D+00
        iside = iside + 1
     else if ( iside == 2 ) then
        v(2) = 0.0D+00
        v(1) = 1.0D+00
        v(4) = 0.0D+00
        iside = iside + 1
     end if
     !
     !  Calculate the error near the boundary layer at 1.
     !
  else if ( mode == 4 ) then

     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) '      X             G(X)          G(X)-F(X):'
     write ( *, '(a)' ) ' '

     x = 0.75D+00

     do i = 1, 9

        ep1 = exp ( s2ovep * ( 1.0D+00 - x ) ) * factor
        ep2 = exp ( s2ovep * ( 1.0D+00 + x ) ) * factor

        solutn = 12.0D+00 / ( 1.0D+00 + ep1 )**2 * ep1 &
             + 12.0D+00 / ( 1.0D+00 + ep2 )**2 * ep2 - 1.0D+00

        value = ppvalu ( break, coef, l, kpm, x, 0 )

        error = solutn - value
        write ( *, '(2x,3g14.6)' ) x, solutn, error
        x = x + 0.03125D+00

     end do

  else

     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) 'DIFEQU - Fatal error!'
     write ( *, '(a)' ) '  Illegal value of MODE:'
     write ( *, '(a,i8)' ) mode
     stop

  end if

  return
end subroutine difequ