A multi-component lumped kinetic model of non-isothermal and non-linear reactive liquid chromatography was formulated and approximated numerically to demonstrate thermal effects on reaction kinetics, adsorption equilibria, and conversion-separation studies in thermally insulated, packed bed, chromatographic reactors. The considered model is constituted of systems of non-linear convection diffusion reaction partial differential equations for mass and energy balances in the bulk phase coupled with differential equations for mass and energy balances in the stationary phase. In this work, a total variation bounded (TVB) Runge-Kutta local-projection discontinuous Galerkin finite element method (DG-FEM) was derived and proposed for the numerical solutions of the model equations. The developed numerical method is robust, explicit, capable of resolving sharp discontinuities and is second-order accurate. System parametric studies treating heterogeneously catalyzed reversible reactions were performed through numerical simulations. The coupling between thermal and concentration fronts, the influence of temperature on reactor efficiency, and the conversion-separation of products are demonstrated through several consistency tests. The results, which authenticate the accuracy of the (DG-FEM) method, will be beneficial for interpreting mass and energy profiles in non-equilibrium and non-isothermal liquid chromatographic reactors and provide deeper insight into the sensitivity of the conversion-separation process.