PROGRAM runmeb USE common703 USE toms703 IMPLICIT NONE INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE REAL (dp) :: y(20) ! THE INCLUSION OF THE COMMON BLOCK /MEBDF2/ ALLOWS THE USER TO ! ACCESS VARIOUS COUNTERS THAT MAY BE OF INTEREST TO HIM. ! THESE VARIABLES ARE EXPLAINED IN THE COMMENTS IN SUBROUTINE MEBDF. INTEGER, PARAMETER :: lout = 6 REAL (dp) :: err1, err2, err3, err4, hstart, h0, tol, x, xend, xout INTEGER :: index, n, nn n = 4 DO nn=3,10 INDEX = 1 xend = 20.0_dp x = 0.0_dp y(1) = 1.0_dp y(2) = 1.0_dp y(3) = 1.0_dp y(4) = 1.0_dp IF (nn == 3 .OR. nn == 4) hstart = 1.0D-4 IF (nn == 5 .OR. nn == 6) hstart = 1.0D-5 IF (nn == 7 .OR. nn == 8) hstart = 1.0D-6 IF (nn == 9 .OR. nn == 10) hstart = 1.0D-7 ! THESE ARE THE INITIAL STEPS CHOSEN BY THE DETEST PACKAGE FOR THIS PROBLEM h0 = hstart tol = 10.0_dp**(-nn) xout = 20.0_dp 10 CALL mebdf(n, x, h0, y, xout, xend, tol, 21, INDEX, lout, fcn) IF ((INDEX /= 0) .AND. (INDEX /= 3)) THEN IF(INDEX == 1) THEN INDEX = 0 GO TO 10 END IF WRITE( 6 ,15) INDEX 15 FORMAT(' ***INTEGRATION HAS FAILED*** WITH INDEX = ', i3) STOP END IF ! HAVE COMPLETED ONE STEP ! HAVE WE FINISHED YET? IF( INDEX == 0) THEN ! THEN WE HAVE EFFECTIVELY HIT TOUT x = xout ELSE INDEX = 3 GO TO 10 END IF WRITE(6, 690) tol 690 FORMAT(' REQUESTED TOLERANCE ', G10.3) WRITE(6,700) nfe, nje 700 FORMAT(' ', i5, ' FUNCTION EVALUATIONS ', & i5, ' JACOBIAN EVALUATIONS') err1 = y(1) - 4.539992969929191D-05 err2 = y(2) - 2.061153036149920D-09 err3 = y(3) - 2.823006338263857D-56 err4 = y(4) - 5.235792540515384D-52 ! THESE VERY SMALL CONSTANTS SHOULD BE SET TO ZERO IF THERE ! ARE LIKELY TO BE DIFFICULTIES DUE TO UNDERFLOW. WRITE(6,710) 710 FORMAT(' ERRORS AT ENDPOINT') WRITE(6,720) err1, err2, err3, err4 720 FORMAT(' ERRS', 5G14.5) END DO STOP END PROGRAM runmeb SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) ydot(1) = -0.5_dp*y(1) ydot(2) = -y(2) ydot(3) = -100.0_dp*y(3) ydot(4) = -90.0_dp*y(4) ! N.B. This application does not use argument t, but for ELF90: IF (t < -9.9E99_dp) STOP RETURN END SUBROUTINE fcn SUBROUTINE pderv(t, y, pw) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:,:) REAL (dp), INTENT(OUT) :: pw(:,:) ! pw(4,4) pw(1:4,1:4) = 0.0_dp pw(1,1) = -0.5_dp pw(2,2) = -1.0_dp pw(3,3) = -100.0_dp pw(4,4) = -90.0_dp ! N.B. This application does not use arguments t & y, but for ELF90: IF (t < -9.9E99_dp) STOP IF (y(1,1) < -9.9E99_dp) STOP RETURN END SUBROUTINE pderv