Look at Neural ODE from Neurips 2018

Link for the lazy: https://arxiv.org/abs/1806.07366

I think this one is a good start:
Differentiation of Blackbox Combinatorial Solvers
https://arxiv.org/abs/1912.02175

I am making a solver in pytorch :p

https://docs.sciml.ai/DiffEqFlux/stable/

Someone has already mentioned Neural Ordinary Differential Equations, which is also the first thing that came to mind. There are also extensions to it, where one can use PDEs(Neural Hamiltonian Flows) or even stochastic DEs(Score-Based Generative Models) in the model. All of them covering different but overlapping use cases.

There are also techniques which use numerical solvers as blackboxes to perform model-order reduction of a complicated system of equations, or identifying slow modes, timescale decomposition, etc.

Might be helpful https://www.pgupta.info/data/talks/huawei23.pdf

Typically, if you're solver is not written in Pytorch/tensorflow itself you can't easily calculate gradients through them as your computational graph doesn't capture the solver. If your soler is also written in the framework and differentiable you might be able to just backpropagate through it though. Otherwise, the Neural ODE paper that was linked here a few times has an adjoint formulation that gives you the gradient wrt to the solver as a solution to another ode, but this is specific to their problem and won't apply to non-differential equations.

I don't know which numerical optimization OP is trying to solve, but one major weakness of this paper is that their method requires two solver calls per instance per epoch... Training time might quickly become intractable.

OP should have a look at other methods that aim to solve their problem efficiently, such as https://arxiv.org/abs/2112.03609 or recently https://arxiv.org/abs/2203.16067