Some more testing is needed, but looks promising so far

Test that all derivatives match with those of finite differentiation to numerical precision (they do!)

Also not clear if this is the best way to include it.  Probably need the std::variant approach here too.