f = @(p, x) ... where x is one row of indeps, and f returns two rows, each
row for one equation (each equation uses all p elements). This code works
OK.
If I wanted to use explicit jacobian to improve precision and calculation
stability, please how should I handle the two equations in f and dfdp?
Sorry, please forget what I said about the reshaping of 'x', I was
confused... I probably shouldn't try to answer help requests while I'm
working on something different. And I probably don't understand what
you mean with 'the 2D case'.
The Jacobian should just compute one row of partial derivatives for
each element of the result of 'f'.