subroutine jacobi(f, fx, x, J)
  implicit none

  integer, parameter           :: m = 2, n = 2
  real (kind = 8), intent(in)  :: fx(m), x(n)
  real (kind = 8), intent(out) :: J(m, n)
  real (kind = 8) :: E(n,n), fxplushe(m), xplushe(n)
  real (kind = 8), parameter   :: h = 1.4901D-08
  integer                      :: i
  logical, parameter           :: printing = .false.

  interface 
     subroutine f(x, fret, printing)
       implicit none
       real(kind = 8), intent(in)  :: x(2)
       real(kind = 8), intent(out) :: fret(2)
       logical, intent(in)         :: printing
     end subroutine f
  end interface

  ! call f(x, fx)

  E(1,1) = 1
  E(1,2) = 0
  E(2,1) = 0
  E(2,2) = 1

  do i = 1, n
     xplushe = x + h*E(:, i)
     call f(xplushe, fxplushe, printing)
     J(:, i) = (fxplushe - fx)/h
  end do
end subroutine jacobi