C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR EXTRAPOLATION CODES SEULEX AT VAN DER POL
C * * * * * * * * * * * * * * * * * * * * * * * * *
clink dr_seulex seulex dc_decsol decsol
clink dr_seulex seulex dc_lapack lapack lapackc
        IMPLICIT REAL*8 (A-H,O-Z)


C --- PARAMETERS FOR SEULEX (FULL JACOBIAN)
        PARAMETER (ND=3,KM=12,KM2=2+KM*(KM+3)/2,NRDENS=ND)
        PARAMETER (LWORK=2*ND*ND+(KM+8)*ND+4*KM+20+KM2*NRDENS)
        PARAMETER (LIWORK=2*ND+KM+20+NRDENS) 
C --- DECLARATIONS
C-------CPU TIME
        real ta(2), t
        real dtime
        external dtime       
        DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),RTOL(ND),ATOL(ND)
        EXTERNAL FVPOL,JVPOL,SOLOUT 
        RPAR=1.0D-6

C --- DIMENSION OF THE SYSTEM
        N=3
C --- PROBLEM IS AUTONOMOUS
        IFCN=1
C --- COMPUTE THE JACOBIAN ANALYTICALLY
        IJAC=0
C --- JACOBIAN IS A FULL MATRIX
        MLJAC=N
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
        IMAS=0
C --- OUTPUT ROUTINE IS USED DURING INTEGRATION
        IOUT=2
C --- INITIAL VALUES
        t = dtime(ta)     ! start timing
        X=0.0D0
        Y(1)=1.0D0
        Y(2)=0.0D0
        Y(3)=0.0D0
C --- ENDPOINT OF INTEGRATION
        XEND=4.0D10
C --- REQUIRED TOLERANCE  
        ITOL=1      
        RTOL(1)=1D-4
	     RTOL(2)=1D-4
	     RTOL(3)=1D-4
        ATOL(1)=1D-10
	     ATOL(2)=1D-15
	     ATOL(3)=1D-3
        
C --- INITIAL STEP SIZE
        H=1.0D-6 
C --- SET DEFAULT VALUES
        DO 10 I=1,20
  10    WORK(I)=0.D0
        DO 12 I=1,20
  12    IWORK(I)=0
        IWORK(6)=N
        DO I=1,NRDENS
          IWORK(20+I)=I
        END DO
C --- CALL OF THE SUBROUTINE SEULEX 
        CALL SEULEX(N,FVPOL,IFCN,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JVPOL,IJAC,MLJAC,MUJAC,
     &                  FVPOL,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
C --- PRINT FINAL SOLUTION
        WRITE (6,99) X,Y(1),Y(2),Y(3)
 99     FORMAT(1X,'X =',E18.10,'    Y =',3E18.10)
C --- PRINT STATISTICS
C        WRITE (6,90) RTOL
C 90     FORMAT('       rtol=',D8.2)
        WRITE (6,91) (IWORK(J),J=14,20)
 91     FORMAT(' fcn=',I5,' jac=',I4,' step=',I4,
     &        ' accpt=',I4,' rejct=',I3,' dec=',I4,
     &        ' sol=',I5)

C---CPU TIME
        t = dtime(ta)     ! stop timing
        print *, 'Elapsed CPU time in seconds =', t
C-----END

        STOP
        END
C
C
        SUBROUTINE SOLOUT (NR,XOLD,X,Y,RC,LRC,IC,LIC,N,
     &                     RPAR,IPAR,IRTRN)
C --- PRINTS SOLUTION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),RC(LRC),IC(LIC)
        COMMON /INTERN/XOUT
        REAL  Y1,Y2,Y3
        IF (X.LE.0.4) THEN
            XOUT=0.4
        END IF
        IF (X.GE.XOUT) THEN
           Y1=CONTEX(1,XOUT,RC,LRC,IC,LIC)
           Y2=CONTEX(2,XOUT,RC,LRC,IC,LIC)
           Y3=CONTEX(3,XOUT,RC,LRC,IC,LIC)
           WRITE (*,99) XOUT,Y1,Y2,Y3
C          WRITE (*,*)  CONTEX(2,XOUT,RC,LRC,IC,LIC),CONTEX(3,XOUT,RC,LRC,IC,LIC)
C          WRITE (*,*)  CONTEX(3,XOUT,RC,LRC,IC,LIC)
            XOUT=XOUT*10.0D0
        END IF
 99     FORMAT(1X,'X =',E18.10,'    Y =',3E18.10)
        RETURN
        END
C
C
        SUBROUTINE FVPOL(N,X,Y,F,RPAR,IPAR)
C --- RIGHT-HAND SIDE OF VAN DER POL'S EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),F(N)
        F(1)=-0.04*Y(1)+1.0D4*Y(2)*Y(3)
        F(3)=3.0D7*Y(2)**2
        F(2)=-F(1)-F(3)
        RETURN
        END 
C
C
        SUBROUTINE JVPOL(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN OF VAN DER POL'S EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),DFY(LDFY,N)
        DFY(1,1)=-0.04
        DFY(1,2)=1.0D4*Y(3)
        DFY(1,3)=1.0D4*Y(2)
        DFY(2,1)=0.04
        DFY(2,2)=-1.0D4*Y(3)-6.0D7*Y(2)
        DFY(2,3)=-1.0D4*Y(2)
        DFY(3,1)=0.0
        DFY(3,2)=6.0D7*Y(2)
        DFY(3,3)=0
        RETURN
        END