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