@Ianast,

Can you please confirm you have checked your Fortran code, especially your function evaluations, using suitable and small unit tests and ensured the oft overlooked items such as precision and arithmetic of the floating-point calculations are all consistent with your needs?

You can review the trivial example below to see the different possibilities.

```
blk1: block
real :: x, dx, dfdx
print *, "Block 1: use default precision and mixed-mode arithmetic"
x = 0.4
dx = sqrt( epsilon(dx) )
dfdx = ( sin(x+dx) - sin(x-dx) )/ (2*dx)
print *, "f:= sin(x); at x=0.4, dfdx = ", dfdx
print *, "analytical solution: ", cos(x)
print "(g0,1x,1pg10.2)", "% error", (dfdx/cos(x) - 1.0)*100.0
end block blk1
print *
blk2: block
integer, parameter :: wp = selected_real_kind( p=12 )
real(kind=wp) :: x, dx, dfdx
print *, "Block 2: use defined precision of ", precision(x), " and consistent arithmetic"
x = 0.4_wp
dx = sqrt( epsilon(dx) )
dfdx = ( sin(x+dx) - sin(x-dx) )/ (2.0_wp*dx)
print *, "f:= sin(x); at x=0.4, dfdx = ", dfdx
print *, "analytical solution: ", cos(x)
print "(g0,1x,1pg10.2)", "% error", (dfdx/cos(x) - 1.0_wp)*100.0_wp
end block blk2
end
```

Upon execution of the program, the expected output is

Block 1: use default precision and mixed-mode arithmetic

f:= sin(x); at x=0.4, dfdx = 0.921042860

analytical solution: 0.921060979

% error -1.97E-03

Block 2: use defined precision of 15 and consistent arithmetic

f:= sin(x); at x=0.4, dfdx = 0.92106099426746368

analytical solution: 0.92106099400288510

% error 2.87E-08

As you know, such small unit tests can help you considerably with debugging the differences.