In this computational paper, we focused on the efficient numerical implementation of semi-implicit methods for models in materials science. In particular, we were interested in a class of nonlinear higher-order parabolic partial differential equations. The Cahn-Hilliard (CH) equation was chosen as a benchmark problem for our proposed methods. We first considered the Cahn-Hilliard equation with a convexity-splitting (CS) approach coupled with a backward Euler approximation of the time derivative and tested the performance against the bi-harmonic-modified (BHM) approach in terms of accuracy, order of convergence, and computation time. Higher-order time-stepping techniques that allow for the methods to increase their accuracy and order of convergence were then introduced. The proposed schemes in this paper were found to be very efficient for 2D computations. Computed dynamics in 2D and 3D are presented to demonstrate the energy-decreasing property and overall performance of the methods for longer simulation runs with a variety of initial conditions. In addition, we also present a simple yet powerful way to accelerate the computations by using MATLAB built-in commands to perform GPU implementations of the schemes. We show that it is possible to accelerate computations for the CH equation in 3D by a factor of 80, provided the hardware is capable enough.