In this thesis, we present numerical solution of semilinear partial differential equations (PDEs) where the linear differential operator is a self-adjoint. A recent spectral method for self-adjoint operators, based on basis recombination, leads to symmetric definite matrices, which have real spectrum. This allows for developing stable time-stepping algorithm to solve the resulting the ordinary differential equations. The linear part of the discretized problem is usually stiff, which constraints the step size for explicit numerical schemes. We therefore use exponential integrators, a well-known time-stepping methods for solving stiff differential equations. We describe three methods namely, eigen-decomposition, contour integral and Carath'{e}odory-Fej'{e}r approximation, for computing the matrix ($\varphi$) functions of the exponential integrators. We perform numerical experiments with some PDEs with different boundary conditions, including time-dependent boundary condition and the numerical results confirm the accuracy of the methods in both space and time.