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,2*NLAGS) :: Z
    DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: YLAG,DYLAG
    INTENT(IN)  :: T,Y,Z
    INTENT(OUT) :: DY
    ! NOTE:
    ! This is a neutral problem. Columns 1...NLAGS of the Z
    ! matrix contain the delayed solution values and columns
    ! NLAGS+1...2*NLAGS contain the delayed derivative values.
     
      YLAG  = Z(:,1:NLAGS)
      DYLAG = Z(:,NLAGS+1:2*NLAGS)

      DY(1) = COS(T)*(1D0 + YLAG(1,1)) + 0.3D0*Y(1)*DYLAG(1,1)  &
        + 0.7D0*SIN(T)*COS( T*SIN(T)**2 ) - SIN( T+T*SIN(T)**2 )
                   
    RETURN
  END SUBROUTINE DDES

  SUBROUTINE HISTORY(T,Y)

    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN,2) :: Y
    INTENT(IN)  :: T
    INTENT(OUT) :: Y
    ! NOTE:
    ! This is a neutral problem.  Return the history
    ! in column 1 of the Y matrix and the derivative
    ! of the history in column 2.

      Y(1,1) = SIN(T)
      Y(1,2) = COS(T)

    RETURN
  END SUBROUTINE HISTORY

  SUBROUTINE BETA(T,Y,BVAL)
    
    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y
    DOUBLE PRECISION, DIMENSION(NLAGS) :: BVAL
    INTENT(IN)  :: T,Y
    INTENT(OUT) :: BVAL

      BVAL(1) = T*Y(1)**2

    RETURN
  END SUBROUTINE BETA

  SUBROUTINE CHECK(T,Y,ERROR)
    
    DOUBLE PRECISION :: T,Y
    DOUBLE PRECISION :: ERROR
    INTENT(IN)  :: T,Y
    INTENT(OUT) :: ERROR

      ERROR = ABS(Y - SIN(T))

    RETURN
  END SUBROUTINE CHECK

END MODULE define_DDEs

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

PROGRAM neutvan 

! Program to solve the last example in the SGT book. This 
! is a neutral problem with a state-dependent vanishing delay.
!
! The DDE is defined in the module define_DDEs.  The problem
! is solved here with DDE_SOLVER and its output written to a
! file.  The auxilary function NEUTVAN.M imports the data into
! Matlab and plots 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 
  ! 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 
  ! 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

  OPTS = DDE_SET(RE=1D-10,AE=1D-10,NEUTRAL=.TRUE.)

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

  ! Was the solver successful?
  IF(SOL%FLAG == 0) THEN

    PRINT *,' '    
    ! Print integration statistics:
    CALL PRINT_STATS(SOL)
 
    ! Write the solution to a file for subsequent 
    ! plotting in Matlab.
    OPEN(UNIT=6, FILE='neutvan.dat')
    DO I = 1,SOL%NPTS
      WRITE(UNIT=6,FMT='(2D12.4)') SOL%T(I),(SOL%Y(I,J),J=1,NEQN)
    ENDDO

    PRINT *,' Normal return from DDE_SOLVER with results'
    PRINT *," written to the file 'neutvan.dat'."
    PRINT *,' '
    PRINT *,' These results can be accessed in Matlab'
    PRINT *,' and plotted by'
    PRINT *,' '
    PRINT *," >> [t,y] = neutvan;"
    PRINT *,' '
    
    MAXERROR = 0D0
    DO I = 1,SOL%NPTS
       CALL CHECK(SOL%T(I),SOL%Y(I,1),ERROR)
       MAXERROR = MAX(ERROR,MAXERROR)
    END DO
    PRINT *,' '
    PRINT *,' Comparison to the analytical solution shows'
    PRINT *,' that at all steps the absolute error is less'
    PRINT *,' than ',MAXERROR,'.'
    PRINT *,' '

  ELSE

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

  ENDIF
  
  STOP
END PROGRAM neutvan