Manual Newton's method faster than fsolve with jacobian for a nonlinear system

I have a system of 100 nonlinear equations with 100 variables.

• I use fsolve (giving it the analytical Jacobian) and it takes on average 0.033 seconds
• I manually do Newton's method with the same tolerance level (1e-10) and it takes 0.009 seconds.

The quality of the results and their sensitivity to starting values are statistically the same. This difference in speed remains similar when I scale up to 1000 equations/variables [0.68 and 0.38 respectively].

Question: Am I misusing fsolve somehow? There's no way a fine-tuned optimizer like fsolve should be worse.

Note: Providing a working example here is too difficult as the function, its inputs, and the Jacobian take over a hundred lines.

I can provide the general syntax I'm using:

options = optimoptions('fsolve','SpecifyObjectiveGradient',true,...
    'MaxIter',10000,'Display','iter','TolFun',1e-6,'TolX',1e-6,'MaxFunEvals',30000);

[x]=fsolve(@(t)function(t,inputs),x0,options);

The manual method is as follows:

x0_new=x0;
dx=1;
while sum(dx.^2)>1e-6
    x0_0=x0_new;
    dx=function_taylor(x0_0,inputs);
    x0_new=x0_0+dx;
end

UPDATE: So I ran the profiler and I think it's the overhead in fsolve that causes the issue. Once I increase the number of equations, the speed differential declines a bit, and this is likely due to the shrinking share of time costs that the fsolve overhead represents.

[Milind Jain]

Can you please send us the code so that we can reproduce the same results at our end. This would help us in better understanding the issue and suggest improvements.