opts = optimset('Jacobian', 'On'); x = fsolve(@fun_ok, ones(3,1), opts) x = fsolve(@fun_wrong, ones(3,1), opts)