MODULE Stiff   ! Defines a stiff system of 2 ODEs
               ! Stiffness increases with x. Jacobian has
               ! eigenvalues -x^4 and -2x^(-0.5)
               ! True solution given
!------------------------------------J R Dormand,   1/1996
! Modified: A Chisalita, 1/2006
! ------------------------------

  INTEGER, PARAMETER :: neq = 2
  INTEGER :: nfunc, njac

CONTAINS
 !

  SUBROUTINE Output(idev, x, y, true_solution)          !  Outputs solution with
    REAL (KIND = 8) :: x, y(:)     !  error using True
logical true_solution
integer idev

if(true_solution) then
write(idev, '(f10.5, 2(1pg19.10), 2(1pg 15.4))')  x, y(:), True(x) -y(:)
else
write(idev, '(f10.5, 2(1pg19.10), 2(1pg 15.4))')  x, y(:)
endif

  END SUBROUTINE Output

!
  FUNCTION Fcn(x, y)             ! ----- Defines derivative
                             ! This stiffens as x increases
    REAL (KIND = 8) :: x, y(:), Fcn(SIZE(y))
      nfunc = nfunc + 1
!      Fcn(1) = 1.0d0/y(1) - x**2 - 2.0d0/x**3
!      Fcn(2) = y(1)/y(2)**2 - 1.0d0/x - 0.5d0/SQRT(x**3)

fcn(1) =y(2)
fcn(2) =-100*y(1) -101*y(2)

  END FUNCTION Fcn

!
  FUNCTION Jac(x, y)             ! - Jacobian matrix of Fcn
    REAL (KIND = 8) :: x, y(:), Jac(neq, neq)
! A simple stiff example (JRD 1/95)
! TEST example:

     njac = njac + 1
!     Jac(1,1) = -1.0d0/y(1)**2
!     Jac(1,2) =  0.0d0
!     Jac(2,1) =  1.0d0/y(2)**2
!     Jac(2,2) = -2.0d0*y(1)/y(2)**3

jac(1,1) =0;	jac(1,2) =1;
jac(2,1) =-100;	jac(2,2) =-101

  END FUNCTION Jac

!
  FUNCTION True(x)                  !  True solution y(x)
    REAL (KIND = 8) :: x, True(neq)
!    True = (/ 1.0d0/x**2,  1.0d0/SQRT(x) /)
True(1) =(100._8*dexp(-x) -1_8*dexp(-100*x) )/99._8
True(2) =100._8/99._8*(-dexp(-x) +dexp(-100*x))

  END FUNCTION True
!
END MODULE Stiff
