Code:

MODULE GLOBVARS
REAL ARR0(3,3) , ARR1(3,3) , DELTIME
END MODULE GLOBVARS
PROGRAM MAIN
INCLUDE 'link_fnl_static.h'
USE GLOBVARS
USE IVPAG_INT
USE UMACH_INT !Routine UMACH sets or retrieves the input, output, or error output device unit numbers.
USE LINRG_INT
IMPLICIT NONE
INTEGER MXPARM, N
PARAMETER (MXPARM=50, N=4)
! SPECIFICATIONS FOR IVPAG PARAMETERS
INTEGER MABSE, MBDF, MSOLVE
PARAMETER (MABSE=1, MBDF=2, MSOLVE=2)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER IDO, ISTEP, NOUT
REAL A(1,1), PARAM(MXPARM), T, TEND, TOL, Y(N)
!
REAL D0(3,3), D1(3,3),DT
!
! SPECIFICATIONS FOR FUNCTIONS
EXTERNAL FCN, FCNJ
!
DATA D0/1.985,0.0,0.0,0.0,0.725,0.0,0.0,0.0,0.7/,
1D1/2.0,0.0,0.0,0.0,0.718,0.0,0.0,0.0,0.695/, DT/0.1/
!
ARR0 = D0
ARR1 = D1
DELTIME = DT
!
!Routine UMACH sets or retrieves the input, output, or error output device unit numbers.
CALL UMACH (2, NOUT) !2 - Retrieves and stores the output unit number in NOUT
! Set initial conditions
T = 0.0
Y(1) = 1.0
Y(2) = 1.0
Y(3) = 1.0
Y(4) = 1.0
! Set TEND
TEND = 10.0
! Set error tolerance
TOL = 0.000001
! Set PARAM to defaults
PARAM = 0.0E0
!
PARAM(10) = MABSE
! Select BDF method
PARAM(12) = MBDF
! Select chord method and a divided difference Jacobian.
PARAM(13) = MSOLVE
! Solve & Print
WRITE (NOUT,99998)
99998 FORMAT (5X, 'Time', 9X, 'Y1', 11X, 'Y2', 11X,'Y3',
111X, 'Y4')
IDO = 1
10 CONTINUE
WRITE(NOUT,'(A)') 'CALLING IVPAG'
CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y, TOL=TOL, PARAM=PARAM)
WRITE(NOUT,'(A)') 'EXITING IVPAG'
IF (T .LE. TEND .AND. IDO .NE. 3) THEN
WRITE (NOUT,'(5F12.6)') T, Y
IF (T .EQ. TEND) THEN
IDO = 3
GO TO 10
END IF
END IF
! Show number of function calls.
WRITE (NOUT,99999) PARAM(35)
99999 FORMAT (4X, 'Number of fcn calls with IVPAG =', I6)
END
SUBROUTINE FCN (N, T, Y, YPRIME)
USE GLOBVARS
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(N), YPRIME(N), A(3,3), B(3,0), INV(3,3), X(3,0),
1FDOT(3,3), FINV(3,3), A1(3,3)
!
WRITE(NOUT,'(A)') 'ENTERING SUBROUTINE FCN'
!
FDOT = (ARR1-ARR0)/DELTIME
A1 = ARR1
!
CALL LINRG (A1, AINV)
FINV = AINV
!
YPRIME(1) = (FDOT(1,1)*FINV(1,1)*Y(1) + FDOT(1,2)*FINV(2,1)*Y(1)
1+FDOT(1,2)*FINV(2,2)*Y(4)+FDOT(1,3)*FINV(3,1)*Y(1))
2+0.5*((3/((Y(2)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))
3+(Y(1)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))+(1/Y(3)))) - Y(1))
!
YPRIME(2) = (FDOT(2,1)*FINV(1,1)*Y(4) + FDOT(2,1)*FINV(1,2)*Y(2)
1+FDOT(2,3)*FINV(3,1)*Y(4) + FDOT(2,3)*FINV(3,2)*Y(2))
2+0.5*((3/((Y(2)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))
3+(Y(1)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))+(1/Y(3)))) - Y(2))
!
YPRIME(3) = (FDOT(3,1)*FINV(1,3)*Y(3) + FDOT(3,2)*FINV(2,3)*Y(3)
1+FDOT(3,3)*FINV(3,3)*Y(3))
2+0.5*(3/((Y(2)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))
3+(Y(1)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4))) + (1/Y(3))) -Y(3))
!
YPRIME(4) = (FDOT(1,1)*FINV(1,1)*Y(4) + FDOT(1,1)*FINV(1,2)*Y(2)
1+FDOT(1,2)*FINV(2,1)*Y(4) + FDOT(1,2)*FINV(2,2)*Y(2)
2+FDOT(1,3)*FINV(3,1)*Y(4) + FDOT(1,3)*FINV(3,2)*Y(2))
RETURN
END
SUBROUTINE FCNJ (N, T, Y, DYPDY)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(N), DYPDY(N,*)
!
RETURN
END

Thank You