MODULE DEFINE_FCN

! Define global variables for the number of differential
! equations, NEQNS, and the number of boundary conditions
! at the left end of the interval, LEFTBC.  The number of
! boundary conditions at the right end is NODE-LEFTBC.
  INTEGER, PARAMETER :: NEQNS=2,LEFTBC=1,RIGHTBC=NEQNS-LEFTBC  

CONTAINS

  SUBROUTINE FSUB(X,Y,F)
	DOUBLE PRECISION :: X,Y(NEQNS),F(NEQNS)

    F = (/ Y(2), -(X*Y(1) - 1D0)*Y(1) /)
    
  END SUBROUTINE FSUB
  
  SUBROUTINE DF(X,Y,PD)
	DOUBLE PRECISION :: X,Y(NEQNS),PD(NEQNS,NEQNS)

    PD = 0D0
    PD(1,2) = 1D0
    PD(2,1) = 1D0 - 2D0*X*Y(1)

  END SUBROUTINE DF    

  SUBROUTINE BCSUB(YA,YB,BCA,BCB)
	DOUBLE PRECISION :: YA(NEQNS),YB(NEQNS),BCA(LEFTBC),BCB(RIGHTBC)

	BCA(1) = YA(2)
	BCB(1) = YB(1) + YB(2)

  END SUBROUTINE BCSUB
  
  SUBROUTINE GUESS_Y(X,Y)
    DOUBLE PRECISION :: X,Y(NEQNS)

    IF (X <= 1.5D0) THEN
      Y(1) = 2D0
      Y(2) = 0D0
    ELSE
      Y(1) = 2D0*EXP(1.5D0 - X)
      Y(2) = - Y(1)
    END IF
    
  END SUBROUTINE GUESS_Y 

END MODULE DEFINE_FCN

!************************************************************

PROGRAM ESI
! U.M. Ascher and R.D. Russell, Reformulation of boundary
! value problems into `standard' form, SIAM Review 23 (1981)
! 238-254 use a BVP from electromagnetic self-interaction
! theory to discuss problems with singular coefficients set
! on an infinite interval.  After some preparation they solve
!
!   u'' + 4*u'/t + (t*u - 1)*u = 0
!
! with u'(0) = 0 and u(L) + u'(L) = 0. The problem is set on
! an infinite interval, so some experimentation is necessary
! to verify that a sufficiently large L has been specified.
! They present results for u(0) when L = 5,8,10,20.  They
! use the initial guess u(t) = 2 for 0 <= t <= 1.5 and 
! u(t) = 2*exp(1.5 - t) for t > 1.5.  We use this guess for
! the first L and thereafter use the solution for one L as
! the guess for the next, extending to the right with the
! value from their guess (which has the right asymptotic
! behavior).

  USE DEFINE_FCN
  USE BVP_M
  
  TYPE(BVP_SOL) :: SOL
  DOUBLE PRECISION :: S(NEQNS,NEQNS),L(4),AR(4),YA(NEQNS),YBNEW(NEQNS)
  INTEGER :: I

! Define the singular term.
  S = 0D0
  S(2,2) = -4D0

  L = (/ 5D0, 8D0, 10D0, 20D0 /)
! Results obtained by Ascher/Russell.  
  AR = (/ 2.03464, 2.11986, 2.11997, 2.11997 /) 
  PRINT *,'  L      BVP_SOLVER      ASCHER/RUSSELL' 
  PRINT *,' '
  DO I = 1,4
    IF (I == 1) THEN
      SOL = BVP_INIT(NEQNS,LEFTBC,(/ 0D0,L(I) /),GUESS_Y)
    ELSE
      ! Use current approximate solution info to extend initial guess.
      ! Extend interval to the right. Left end is fixed at 0.
      CALL BVP_EVAL(SOL,0D0,YA)  ! Get solution at 0.0D0.    
      CALL GUESS_Y(L(I),YBNEW)     ! Get solution at L(I).
      SOL = BVP_EXTEND(SOL,0D0,L(I),YA,YBNEW)
    END IF

  ! The code can solve this problem with finite difference Jacobians,
  ! but it is easy to provide the analytical Jacobian and doing so
  ! helps confirm that the code handles this option properly when 
  ! the BVP is singular.
    SOL = BVP_SOLVER(SOL,FSUB,BCSUB,SINGULARTERM=S,DFDY=DF)
 
    CALL BVP_EVAL(SOL,0.0D0,YA)
    WRITE(*,FMT="(F5.0,F14.5,F17.5)") L(I),YA(1),AR(I)
  END DO

  CALL BVP_TERMINATE(SOL)

END PROGRAM ESI
