If you’re using 32 bit floats, that’s way too small of an epsilon. You may want to see this introduction to automatic differentiation which describes why in detail:
So notice that:
julia> eps(Float32)
1.1920929f-7
which means that your derivative term, when 1e-6 is effectively “holding” the value in its 7th decimal place onwards, but you only have 7 decimal places, and thus the accuracy of your derivative will not be more than one digit. So from quick pen and paper, it’s not surprising it’s different in the second decimal place. You’ll do a bit better then to use half of the digits for the differentiation, or in other words, use epsilon as:
julia> sqrt(eps(Float32))
0.00034526698f0
If you really need to use Float32, you may want to use central differencing to improve this a bit (and change epsilon to cbrt(eps(Float32)))