Partial Differential Equations are ubiquitous in scientific computing. While numerical solvers allow for efficiently computing solutions to these PDEs, these solvers can often encounter a variety of issues that stem from both the underlying physical model (e.g., the possibility of shock wave formation) as well as issues caused by the approximation (e.g., introduction of spurious oscillations). These issues can cause the solver’s program execution to either crash (due to overflow from shock formation) or return results with unacceptable levels of inaccuracy (due to spurious oscillations or dissipations). Hence there exists a critical need to apply program analysis to PDE solvers to soundly certify such errors do not arise. Complicating this task is the fact that many PDEs are nonlinear, which in the case of hyperbolic PDEs, is the root cause for the aforementioned pathologies of shock formation, spurious oscillations, and others. As a solution, we develop Phocus, which is the first technique to statically analyze hyperbolic PDE solvers and the first technique to certify precise bounds on nonlinear PDE solutions and other key quantities such as the CFL condition and a solution’s total variation. Hence Phocus can certify the absence of shock formation in the PDE, the stability of the solver, and bounds on the amount of spurious numerical effects. To do so, Phocus uses a novel optimization-based procedure to synthesize precise abstract transformers to enable effective numerical abstract interpretation of finite difference solvers for hyperbolic PDEs. To evaluate Phocus, we develop a novel set of PDE benchmarks and perform an extensive experimental evaluation which demonstrates Phocus’s significant precision and scalability (up to hundreds of time steps)