# Help needed to calculate Jacobian by forward difference method

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
``````

I don’t know what’s the range of your variable, but this line seems wrong: I think you want to have an absolute threshold only if the relative value falls BELOW what can be represented on your machine:

``````if (h==0) h = EPS
``````

In other words,

``````h = max(EPS, EPS*abs(temp))
``````

Another thing I’ve noticed is that you never initialize the unperturbed solution, `fvec`: you should place an initial call at the beginning of `FdJac`:

``````
! Initialize unperturbed vector
call Myfun(n,x,fvec)

do j=1,n
...
end do
``````

its if(h .eq. 0) h = EPS. In call Myfun how can i include do loop if I have multivariable equation. can you brief me.

when include call Myfun(n,x,fvec) at the bigining of FdJac before
do j ) 1,n

the programm is not even compiling

Why not use a library?:

1 Like