1C * * * * * * * * * * * * * * * * * * * * * * * * *
2C --- DRIVER FOR RADAU5 AT VAN DER POL'S EQUATION
3C * * * * * * * * * * * * * * * * * * * * * * * * *
4C link dr1_radau5 radau5 decsol dc_decsol   or
5C link dr1_radau5 radau5 lapack lapackc dc_lapack
6        IMPLICIT REAL*8 (A-H,O-Z)
7C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN)
8        PARAMETER (ND=2,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20)
9        DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK)
10        EXTERNAL FVPOL,JVPOL,SOLOUT
11C --- PARAMETER IN THE DIFFERENTIAL EQUATION
12        RPAR=1.0D-6
13C --- DIMENSION OF THE SYSTEM
14        N=2
15C --- COMPUTE THE JACOBIAN ANALYTICALLY
16        IJAC=1
17C --- JACOBIAN IS A FULL MATRIX
18        MLJAC=N
19C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
20        IMAS=0
21C --- OUTPUT ROUTINE IS USED DURING INTEGRATION
22        IOUT=1
23C --- INITIAL VALUES
24        X=0.0D0
25        Y(1)=2.0D0
26        Y(2)=-0.66D0
27C --- ENDPOINT OF INTEGRATION
28        XEND=2.0D0
29C --- REQUIRED TOLERANCE
30        RTOL=1.0D-4
31        ATOL=1.0D0*RTOL
32        ITOL=0
33C --- INITIAL STEP SIZE
34        H=1.0D-6
35C --- SET DEFAULT VALUES
36        DO I=1,20
37           IWORK(I)=0
38           WORK(I)=0.D0
39        END DO
40C --- CALL OF THE SUBROUTINE RADAU5
41        CALL RADAU5(N,FVPOL,X,Y,XEND,H,
42     &                  RTOL,ATOL,ITOL,
43     &                  JVPOL,IJAC,MLJAC,MUJAC,
44     &                  FVPOL,IMAS,MLMAS,MUMAS,
45     &                  SOLOUT,IOUT,
46     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
47C --- PRINT FINAL SOLUTION
48        WRITE (6,99) X,Y(1),Y(2)
49 99     FORMAT(1X,'X =',F5.2,'    Y =',2E18.10)
50C --- PRINT STATISTICS
51        WRITE (6,90) RTOL
52 90     FORMAT('       rtol=',D8.2)
53        WRITE (6,91) (IWORK(J),J=14,20)
54 91     FORMAT(' fcn=',I5,' jac=',I4,' step=',I4,' accpt=',I4,
55     &        ' rejct=',I3,' dec=',I4,' sol=',I5)
56        STOP
57        END
58C
59C
60        SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN)
61C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5"
62        IMPLICIT REAL*8 (A-H,O-Z)
63        DIMENSION Y(N),CONT(LRC)
64        COMMON /INTERN/XOUT
65        IF (NR.EQ.1) THEN
66           WRITE (6,99) X,Y(1),Y(2),NR-1
67           XOUT=0.2D0
68        ELSE
69 10        CONTINUE
70           IF (X.GE.XOUT) THEN
71C --- CONTINUOUS OUTPUT FOR RADAU5
72              WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC),
73     &                     CONTR5(2,XOUT,CONT,LRC),NR-1
74              XOUT=XOUT+0.2D0
75              GOTO 10
76           END IF
77        END IF
78 99     FORMAT(1X,'X =',F5.2,'    Y =',2E18.10,'    NSTEP =',I4)
79        RETURN
80        END
81C
82C
83        SUBROUTINE FVPOL(N,X,Y,F,RPAR,IPAR)
84C --- RIGHT-HAND SIDE OF VAN DER POL'S EQUATION
85        IMPLICIT REAL*8 (A-H,O-Z)
86        DIMENSION Y(N),F(N)
87        F(1)=Y(2)
88        F(2)=((1-Y(1)**2)*Y(2)-Y(1))/RPAR
89        RETURN
90        END
91C
92C
93        SUBROUTINE JVPOL(N,X,Y,DFY,LDFY,RPAR,IPAR)
94C --- JACOBIAN OF VAN DER POL'S EQUATION
95        IMPLICIT REAL*8 (A-H,O-Z)
96        DIMENSION Y(N),DFY(LDFY,N)
97        DFY(1,1)=0.0D0
98        DFY(1,2)=1.0D0
99        DFY(2,1)=(-2.0D0*Y(1)*Y(2)-1.0D0)/RPAR
100        DFY(2,2)=(1.0D0-Y(1)**2)/RPAR
101        RETURN
102        END
103