It's usually pretty hard to get statistical information about a solution for general problems. I imagine most people just measure statistical variation of the solution by running repeated simulations.
The methods you're talking about usually assume that all of the following are true of the problem
- It is unconstrained and differentiable
- It uses a least squares objective function
- The residuals are Gaussian distributed.
If that is your situation, you could compute the Jacobian at the solution yourself. Or, once you've found the solution, you can feed it to LSQCURVEFIT as an initial point and use its Jacobian output in the usual way. LSQCURVEFIT should stop in 1 iteration, since your solution is already optimal.