jac Subroutine

   pure subroutine jac(t, y, ydot, pd, cj, rpar, ipar)
   !! Jacobian routine.

      real(rk), intent(in) :: t
      real(rk), intent(in) :: y(neq)
      real(rk), intent(in) :: ydot(neq)
      real(rk), intent(out) :: pd(nrowpd, neq)
      real(rk), intent(in) :: cj
      real(rk), intent(in):: rpar
      integer, intent(in) :: ipar

      ! First define the Jacobian matrix for the right-hand side of the ODE:
      ! Y' = F(T,Y), i.e. dF/dY.
      pd(1, 1) = zero
      pd(1, 2) = one
      pd(2, 1) = -200*y(1)*y(2) - one
      pd(2, 2) = 100*(one - y(1)**2)

      ! Next update the Jacobian with the right-hand side to form the DAE Jacobian:
      ! cj*dG/dY' + dG/dY = cj*I - dF/dY.
      pd(1, 1) = cj - pd(1, 1)
      pd(1, 2) = -pd(1, 2)
      pd(2, 1) = -pd(2, 1)
      pd(2, 2) = cj - pd(2, 2)

   end subroutine jac