MODULE define_DDEs

  IMPLICIT NONE
  
  INTEGER, PARAMETER :: NEQN=1,NLAGS=1

CONTAINS

  SUBROUTINE DDES(T,Y,Z,DY)

    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY
    DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: Z

      DY(1) = Y(1) * Z(1,1) / T

    RETURN
  END SUBROUTINE DDES

  SUBROUTINE BETA(T,Y,BVAL)

    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y
    DOUBLE PRECISION, DIMENSION(NLAGS) :: BVAL

      BVAL(1) = LOG(Y(1))

    RETURN
  END SUBROUTINE BETA

  SUBROUTINE CHECK(T,Y,ERROR)

    DOUBLE PRECISION :: T,Y
    DOUBLE PRECISION :: E1, E2, E3, EXACT, ERROR

      E1 = EXP(1D0)
      E2 = EXP(2D0)
      E3 = EXP(3D0 - EXP(1D0 - E1))
      IF (T <= E1) THEN
         EXACT = T
      ELSE
         IF (T <= E2) THEN
            EXACT = EXP(T/E1)
         ELSE
           IF (T <= E3) THEN
              EXACT = (E1/(3D0 - LOG(T)))**E1
           ELSE
              PRINT *, ' Solution not known for T > E3.'
              EXACT = Y
           ENDIF
         ENDIF
      ENDIF
      ERROR = ABS((Y - EXACT)/EXACT)

    RETURN
  END SUBROUTINE CHECK
END MODULE define_DDEs

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

PROGRAM exsd1 

! For t > 1, solve the DDE  y'(t) = y(t)*y(ln(y(t)))/t for t >= 1,
!                           y(t) = 1 for t <= 1.       
!
! The DDE is defined in the module define_DDEs. The problem is
! solved here with DKLAG6 and its output written to a file.
! Using the auxilary function getplotdata, it is then easy to
! import the data into Matlab and plot it.

  USE define_DDEs
  USE DDE_SOLVER_M

  IMPLICIT NONE

  ! The quantities
  !
  !   NEQN  = number of equations
  !   NLAGS = number of delays
  !
  ! are defined in the module define_DDEs as PARAMETERs so
  ! that they can be used for dimensioning arrays here.
  ! They are passed to the solver in the array NVAR.

  INTEGER, DIMENSION(2) :: NVAR = (/NEQN,NLAGS/)

  TYPE(DDE_SOL) :: SOL ! Solution in Matlab style. 
  ! The fields of SOL are expressed in terms of the
  ! number of differential equations, NEQN, and the 
  ! number of output points, NPTS:

  !   SOL%NPTS         -- NPTS,number of output points.
  !
  !   SOL%T(NPTS)      -- values of independent variable, T.  
  !
  !   SOL%Y(NPTS,NEQN) -- values of dependent variable, Y,
  !                       corresponding to values of SOL%T. 

  TYPE(DDE_OPTS) :: OPTS

  ! Local variables:
  INTEGER :: I,J
  DOUBLE PRECISION :: ERROR,MAXERROR,TVAL,YVAL

  CHARACTER(7+6*NEQN) :: EXPORT
  ! FNAME is the name of the output file.
  CHARACTER(10) :: FNAME='export.dat'  

  OPTS = DDE_SET(RE=1D-5,AE=1D-5)

  SOL = DDE_SOLVER(NVAR,DDES,BETA,HISTORY=(/ 1D0 /),&
                  TSPAN=(/ 1D0, 10D0 /),OPTIONS=OPTS)

  ! Was the solver successful?
  IF (SOL%FLAG == 0) THEN
 
     ! Write the solution to a file for subsequent plotting
     ! in Matlab.
     OPEN(UNIT=6, FILE=FNAME)
     ! Set up format for exporting data.
     EXPORT = "(D12.4"//REPEAT(",D12.4",NEQN)//")"
     DO I = 1,SOL%NPTS
        WRITE(UNIT=6,FMT=EXPORT) SOL%T(I),(SOL%Y(I,J),J=1,NEQN)
     ENDDO
      
     PRINT *,' Normal return from DDE_SOLVER with results'
     PRINT *," written to the file '",FNAME,"'."
     PRINT *,' '
     PRINT *,' These results can be accessed in Matlab'
     PRINT *,' from the folder holding this file by'
     PRINT *,' '
     PRINT *," >> [t,y] = getplotdata('",FNAME,"');"
     PRINT *,' '
     PRINT *,' and then plotted with, e.g., '
     PRINT *,' '
     PRINT *,' >> plot(t,y)'
     PRINT *,' '

     MAXERROR = 0D0
     DO I = 1,SOL%NPTS
        TVAL = SOL%T(I)
        YVAL = SOL%Y(I,1)
        CALL CHECK(TVAL,YVAL,ERROR)
        MAXERROR = MAX(ERROR,MAXERROR)
     ENDDO
     PRINT *,' '
     PRINT *,' Comparison to the analytical solution shows that'
     PRINT *,' at all steps the error is less than ',MAXERROR
     PRINT *,' '

  ELSE

     PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ',&
     SOL%FLAG

  ENDIF
  
  STOP
END PROGRAM exsd1