We describe a new numerical method for elliptic and parabolic partial differential equations, with a wide range of application, from heat conduction, fluid dynamics, acoustic propagation, electrostatics, to the option pricing problem in mathematical finance. Unlike approximations of the solution through discretization by the existing methods (such as finite difference, finite element, etc), our method is designed to obtain and optimize hard upper and lower bounding functions of the solution over the domain by mathematical programming, that form a 100% confidence interval within which the solution is guaranteed to exist. An optimization problem is formulated based upon the probabilistic representation of the solution and can be solved by semidefinite programming with the help of sum-of square relaxations. We discuss convergence of hard bounds in smooth problem settings. Numerical results are presented to support our theoretical developments and to illustrate the effectiveness of our methodology and techniques.