% Create some 'magic' functions and data that reproduce the problem. b_pos =[-19.18376 -12.56013;-16.56066 -8.14734]; z =[-68.625;-78.000]; p=[1.8116 -55.0377];s=2.6883; d=@(pos) sqrt( sum((b_pos-pos).^2,2)); Zp=@(pos) permute(sum(permute(d(pos),[1 3 2]).^[columns(p)-1:-1:0].*permute(p,[3 2 1]),2),[1 3 2]); l=@(pos) -sum(log(stdnormal_pdf( (z-Zp(pos)) ./ s ) ./ s)); %Minimize l octave:8> [x, fval, info, output, grad, hess]=fminunc(l,[0 0]);info info = 3 octave:9> grad grad = -0.039276 -0.018594 octave:10> %Note that __fdjac__ is an internal function used to compute the gradient. You must "expose" it somehow octave:10> __fdjac__(l,x,l(x),[1 1],1) ans = -0.032104 -0.022847 octave:11> %As you can see the gradients are different. Now show the hessian octave:11> hess hess = 228.309 -130.248 -130.248 75.519 octave:12> octave:12> %Now try the fix octave:12> [x, fval, info, output, grad, hess]=fminunc2(l,[0 0]);info info = 3 octave:13> grad grad = -0.032104 -0.022847 octave:14> hess hess = 224.124 -133.696 -133.696 80.927