To make sure I understand you correctly – you mean comparing my analytic Jacobian to an AD Jacobian?
If so, comparing maximum(abs.(jac_AD-jac_analytic)) for a random y_in gives me an error of about 1e-7. If I set f=1, this error vanishes, so it’s almost certainly coming from the evaluation of the max function (and either the smoothed or non-smoothed version give similar errors).
However, the MethodError from above persists whether I use max as-is or the smoothed version (both included in the above code). I also tried to smooth the max with log1p(exp(x/epsilon))*epsilon, but x/epsilon can be very large and positive, so this didn’t work as well.
Edit: I just tried some other random y_in vectors and the difference gets up to 1e-3 – still stemming from the max function