C  Test example for the call of subroutine METAN1
C
C  For stiff systems of ordinary differential equations
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C  Number of Equations
      PARAMETER (N      = 4          ,
     2           LOUT   = 6          )
      DIMENSION Y(N)
      COMMON /STATP/ NFCN, NSTEP, NACCPT, NREJCT, NDEC, NSOL
      EXTERNAL FCN
C
C  Print parameter KFLAG
C  KFLAG=0 : No output
C  KFLAG=1 : Integration monitor
C  KFLAG=2 : Intermediate solution points T,Y(1),...
C  KFLAG=3 : Integration monitor and intermediate solution points
      KFLAG=1
C  Starting point of integration
      T=0.D0
C  Initial values Y1(T),...
      DO 1 I=1,4
1     Y(I)=1.D0
C  Final point of integration
      TEND=20.D0
C  Desired accuracy
      TOL=1.D-8
C  Maximum permitted stepsize
      HMAX=TEND-T
C  Initial stepsize guess
C  If H=0, METAN1 internally generates an initial stepsize guess
      H=1.D-3
C
      CALL METAN1 (N,FCN,T,Y,TEND,TOL,HMAX,H,KFLAG)
C
      WRITE(6,101) T,(Y(I),I=1,N)
101   FORMAT(' Solution at T=  ',E25.16,/,10X,4E25.16)
      WRITE (LOUT, '(/,1X,A12, 4(1X, A8), 1X,A6)')
     2      'Tolerance', 'NFCN', 'NDEC',
     3      'NSOL', 'NSTEP', 'KFLAG'
      WRITE (LOUT, '(/,1X,D12.4, 4(1X, I8), 1X,I6)')
     2      TOL, NFCN, NDEC, NSOL, NSTEP, KFLAG
      STOP
      END
      SUBROUTINE FCN (N,T,Y,DY,IFAIL)
C  Example A1 from Stiff Detest(Due to Enright Et Al.)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER IFAIL, N
      DIMENSION Y(N),DY(N)
C
      DY(1)=-0.5D0*Y(1)
      DY(2)=-Y(2)
      DY(3)=-100.D0*Y(3)
      DY(4)=-90.D0*Y(4)
C
      IFAIL=0
      RETURN
      END