In this paper, we use a time splitting method with higher-order accuracy for the solutions (in space variables) of a class of two-dimensional semi-linear parabolic equations. Galerkin-Chebyshev pseudo spectral method is used for discretization of the spatial derivatives, and implicit Euler method is used for temporal discretization. In addition, we use this novel method to solve the well-known semi-linear Poisson-Boltzmann (PB) model equation and obtain solutions with higher-order accuracy. Furthermore, we compare the results obtained by our method for the semi-linear parabolic equation with the available analytical results in the literature for some special cases, and found excellent agreement. Furthermore, our new technique is also applicable for three-dimensional problems.