So I got this working on R2021a, which I'll reproduce here since I've seen other people wonder how to make sfdnls or finitedifferences work and it is undocumented and I don't think anyone has posted an answer previously:
function J = jacobian_matlab(x0, tf)
finDiffOpts.DiffMinChange = 0;
finDiffOpts.DiffMaxChange = inf;
finDiffOpts.TypicalX = ones(length(x0),1);
finDiffOpts.FinDiffType = 'forward';
finDiffOpts.GradObj = 'off';
finDiffOpts.GradConstr = 'off';
finDiffOpts.UseParallel = false;
finDiffOpts.FinDiffRelStep = sqrt(eps);
finDiffOpts.Hessian = [];
sizes.nVar = length(x0);
sizes.mNonlinIneq = 0;
sizes.mNonlinEq = 0;
sizes.xShape = size(x0);
finDiffFlags.fwdFinDiff = true;
finDiffFlags.scaleObjConstr = false;
finDiffFlags.chkComplexObj = false;
finDiffFlags.isGrad = false;
finDiffFlags.hasLBs = false(sizes.nVar,1);
finDiffFlags.hasUBs = false(sizes.nVar,1);
finDiffFlags.chkFunEval = false;
finDiffFlags.sparseFinDiff = false;
funfcn = optimfcnchk(@(x) propagate(x, tf),'fminunc',0,false,false,false)
f = feval(funfcn{3},x0)
J = sfdnls(x0,f,[],[],[],funfcn,[],[],finDiffOpts,sizes,finDiffFlags)
end
The results that I got were that it is about 10x faster than jacobianest, and for my purposes has around similar loss of accuracy at the same point (although the determinant drifts in the opposite direction). I found that at 1e-9 tolerances that numerically integrating the jacobian with ode45 produces much more accurate results in about the same runtime (with higher tolerances to ode45 though it starts to take much more time).
The "propagate" function there is just a function handle to my wrapper around an ode45 problem which is the thing that I'm trying to create a jacobian of.
Probably makes more sense in retrospect to test a full problem rather than to microbenchmark the numerical jacobians, but it looks like numerically integrating the jacobian with ode45 may be a more accurate/stable and equally performant approach for me at least to try.
[ And note that if anyone else tries to get this to work in the future and this doesn't work on future versions of Matlab that you're probably on your own, I'm not providing support for figuring this out, you'll need to consult the sources of stuff like fminunc to see how the options/flags/sizes and the construction of the interior wrapper function happens in your version ]
