Systems of ordinary differential equations (ODEs) and partial differential equations (PDEs) are extensively used in many fields of science, including physics, biochemistry, nonlinear control, and dynamical systems. On the one hand, analytical methods for solving systems of ODEs/PDEs mostly remain an art and are largely insufficient for complex systems. On the other hand, numerical approximation methods do not yield a viable analytical form of the solution that is often required for downstream tasks. In this paper, we present an approximate approach for solving systems of ODEs/PDEs analytically using solvers like Gurobi developed in Operations Research (OR). Our main idea is to represent entire functions as Bézier curves/surfaces with to-be-determined control points. The ODEs/PDEs as well as their boundary conditions can then be reformulated as constraints on these control points. In many cases, this reformulation yields quadratic programming problems (QPPs) that can be solved in polynomial time. It also allows us to reason about inequalities. We demonstrate the success of our approach on several interesting classes of ODEs/PDEs.