The energy balance ordinary differential equations (ODEs) model of climate change is extended to the partial differential equations (PDEs) model with convections and q-diffusions. Instead of integer order second-order partial derivatives, partial q-derivatives are considered. The local stability analysis of the ODEs model is established using the Routh-Hurwitz criterion. A numerical scheme is constructed, which is explicit and second-order in time. For spatial derivatives, second-order central difference formulas are employed. The stability condition of the numerical scheme for the system of convection q-diffusion equations is found. Both types of ODEs and PDEs models are solved with the constructed scheme. A comparison of the constructed scheme with the existing first-order scheme is also made. The graphical results show that global mean surface and ocean temperatures escalate by varying the heat source parameter. Additionally, these newly established techniques demonstrate predictability.