In this paper we develop a method to generate the Lyapunov function for stability analysis for chemical reaction networks. Based on the Chemical Master Equation, we derive the Lyapunov Function partial differential equations (PDEs), whose solution approximates the scaling non-equilibrium potential and serves as the candidate Lyapunov function for the given network. We further prove that for any chemical reaction network the solution (if exists) of the PDEs is dissipative. Moreover, the proposed method of Lyapunov Function PDEs is qualified for analyzing the asymptotic stability of complex balanced networks, all networks with 1-dimensional stoichiometric subspace and some special networks with more than 2-dimensional stoichiometric subspace if some moderate conditions are added. Several examples are presented to illustrate the efficiency of the method.1. Introduction. Chemical reactions networks (CRNs) arise abundantly in the fields including chemistry, systems biology, process industry, and even those seemingly irrelevant to chemistry such as mechanics and ecology. The dynamics of a CRN often appears to be extremely complex due to chemical interactions of cellular processes but sometimes still exhibiting certain regular behaviors like period solutions and stable-fixed-points. As a special subclass, mass-action CRNs (CRNs assigned mass action kinetics, often named mass action systems) have the dynamics of the concentrations of the various species captured by polynomial ordinary differential equations (ODEs), and have received much attention since the pioneering work [10-12,18] emerged. A major concern over this class of systems is to understand the relations between network structures and/or parameters and dynamical properties [6,7,28], especially in characterizing the stability (in the sense of Lyapunov) property [1,2,18,24,27,29]. Following this line of study, we also focus on capturing stability of equilibria in mass action systems (MASs) in the current work. Naturally, Lyapunov functions are desirable objects to prove stability of equilibria. In the field of CRNs, one important example is the pseudo-Helmholtz free energy function, proposed by Horn and Jackson [18]. This Lyapunov function can be derived from the microscopic level using potential theory [2]. Here, we build further on this, and provide a general theory that bridges between the microscopic and the macroscopic level, thermodynamics and potential theory. Based on it, it is possible to derive or find a Lyapunov function for any MAS.The early work on stability analysis mainly focused on exploring causal association from the network topology to the distribution of equilibria, and further to stability of equilibria. Thereinto, the weakly reversible structure, a requirement of complex balanced MAS, is the * most active one. Horn et al. [18] proved the well-known Deficiency Zero Theorem that states a weakly reversible deficiency zero MAS to be complex balanced and to have only one equilibrium in each positive stoichiometric compatibility class. Moreover, ...