I am trying to calculate the Jacobian of the two non linear equations using forward diff method. Here in this code I written a subroutine for the Jacobian. and function for the linear equations.

and main program to write result. when I execute the program I am getting wrong values for the function and Jacobian. I request to explain the code why its happening.

will attach the code below

real function Myfun(x, n, f)

```
implicit none
integer :: n
real :: x(n)
real :: f(n)
f(1) = x(1)**2 + x(2) - 11
f(2) = x(1) + x(2)**2 - 7
```

end function Myfun

program mainMNewton

```
implicit none
real, external :: Myfun
integer :: n
parameter(n=2)
integer :: np
parameter (np=4)
real, dimension(np) :: df
real, dimension(n) :: x
x = reshape((/3.2, -1.7/), shape(x))
call FdJac(n, x, Myfun, np, df)
write(*,*) "Jac:", df
```

end program

```
subroutine FdJac(n, x, fvec, np, df)
integer :: n, np ,Nmax
real :: df(np,np),fvec(n), x(n), EPS
parameter (Nmax = 40, EPS = 1.e-4)
integer :: i, j
real :: h, temp, f(Nmax)
do j = 1, n
temp = x(j)
h = EPS * abs(temp)
if (h .eq. 0) h = EPS
x(j) = temp + h
h = x(j) - temp
call Myfun(n, x, f) ! call for the function vector
x(j) = temp
do i = 1, n
df(i,j) = (f(i)-fvec(i)) / h
write(*,*) "f= ", f(i), "fvec", fvec(i)
end do
end do
return
end subroutine
```