In some previous works, two of the authors have introduced a strategy to develop high-order numerical methods for systems of balance laws that preserve all the stationary solutions of the system. The key ingredient of these methods is a well-balanced reconstruction operator. A strategy has been also introduced to modify any standard reconstruction operator like MUSCL, ENO, CWENO, etc. in order to be well-balanced. This strategy involves a non-linear problem at every cell at every time step that consists in finding the stationary solution whose average is the given cell value. So far this strategy has been only applied to systems whose stationary solution are known either in explicit or implicit form. The goal of this paper is to present a general implementation of this technique that can be applied to any system of balance laws. To do this, the nonlinear problems to be solved in the reconstruction procedure are interpreted as control problems: they consist in finding a solution of an ODE system whose average at the computation interval is given. These problems are written in functional form and the gradient of the functional is computed on the basis of the adjoint problem. Newton’s method is applied then to solve the problems. Special care is put to analyze the effects of computing the averages and the source terms using quadrature formulas. To test their efficiency and well-balancedness, the methods are applied to a number of systems of balance laws, ranging from easy academic systems consisting of Burgers equation with some nonlinear source terms to the shallow water equations or Euler equations of gas dynamics with gravity effects.
