Solved by tom in post #2 
              
              
                So there aren’t efficient ways to integrate similar to autodiff, so basically, you need (more refined) formulas like yours. 
I haven’t seen implementations of the advanced quadrature rules in PyTorch. 
It’s not quite what you want, but the Neural ODE code has ODE solvers which would give you some ad…