      USE DAESL_INT
      USE MP_TYPES
      IMPLICIT NONE 
  
!                             NEQ = Number of equations
      INTEGER, PARAMETER :: NEQ=101
      REAL T, Y(NEQ), YPRIME(NEQ), TEND, X, TRUE, HX, ERRMAX
      INTEGER NOUT, IDO, I, NSTEPS
      REAL, TARGET :: RPARAM(1)
      INTEGER, TARGET :: IPARAM(1)
      TYPE (S_FCN_DATA) USER_DATA
      EXTERNAL GCN
      
!                             Pass NEQ, HX to GCN 
      HX = 1.0 / (NEQ-1)
      IPARAM(1) = NEQ 
      RPARAM(1) = HX 
      USER_DATA%RDATA=>RPARAM
      USER_DATA%IDATA=>IPARAM
!                             Initial values for y, initial guesses for y'
      DO I = 1, NEQ
        X = (I-1) * HX 
        Y(I) = 1 + X 
      END DO
      YPRIME = 0.0
      
      NSTEPS = 10
!                             Always set IDO=1 on first call
      IDO = 1
      DO I = 1, NSTEPS
!                             Output solution at T=0.1,0.2,...,1.0
        T = 0.1 * (I-1)
        TEND = 0.1 * I 
        
!                             Set IDO = 3 on last call
        IF (I == NSTEPS) IDO = 3
        
!                             User-supplied Jacobian matrix (IUJAC=1)
!                             Banded Jacobian (MATSTR=1) 
        CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, IYPR=1, IUJAC=1, &
               MATSTR=1, ML=1, MU=1, RTOL=1.0E-4, FCN_DATA=USER_DATA)
      END DO 
      
      ERRMAX = 0.0
      DO I = 1, NEQ
         X = (I-1) * HX 
         TRUE = (1+X) * EXP(T) 
         ERRMAX = MAX(ERRMAX, ABS(Y(I) - TRUE))
      END DO
      
      CALL UMACH(2, NOUT)
      WRITE (NOUT, *) ' Max Error at T=1 is ', ERRMAX
      END
  
  
      SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)
      USE MP_TYPES
      IMPLICIT NONE
      
      REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), HX
      INTEGER IRES, LDD, I, J, NEQ, MU
      TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA
      
      NEQ = FCN_DATA%IDATA(1)
      HX = FCN_DATA%RDATA(1)
      MU = 1
      
      SELECT CASE (IRES)
      
!                             F_I defined here
      CASE(1)
        DELTA(1) = Y(1) - EXP(T)
        DO I = 2, NEQ-1
           DELTA(I) = -YPRIME(I) + (Y(I+1) - 2.0 * Y(I) + Y(I-1)) &
                      / HX**2 + Y(I)  
        END DO
        DELTA(NEQ) = Y(NEQ) - 2.0 * EXP(T)
        
!                             D(I-J+MU+1,J) = D(F_I)/D(Y_J)
!                             in band storage mode 
      CASE(2)
        D(MU+1,1) = 1.0
        DO I = 2, NEQ-1
           J = I-1
           D(I-J+MU+1, J) = 1.0 / HX**2
           J = I
           D(I-J+MU+1, J) = -2.0 / HX**2 + 1.0
           J = I+1
           D(I-J+MU+1, J) = 1.0 / HX**2
        END DO
        D(MU+1, NEQ) = 1.0
        
!                             D(I-J+MU+1,J) = D(F_I)/D(YPRIME_J)
      CASE(3)
        DO I = 2, NEQ-1
           D(MU+1, I) = -1.0
        END DO

      END SELECT
      END
