As currently programmed, fsolve() does not do numerical derivatives for systems of equations. The reason for this is that subs() is used instead of eval() when evaluating derivatives. [Note: jacob is the symbolically defined jacobian of the system of equations and lsub is the sequence of appropriate numerical substitutions (not the "list" of substitutions the mnemonic might suggest).] The original statement (in `fsolve/sysnewton`) is:
The functions evalf() and subs() apparently do no know how to work together to produce numerical derivatives. However, the following statement does work.
This combination of evalf() and eval() recognizes the need for a numerical derivative and sets up a call to fdiff(). See the help page for fdiff().