{In this study, we present a second-order time-accurate unconditionally stable numerical method for a gradient flow for the Modica–Mortola functional with an equispaced multiple well potential. The proposed second-order time-accurate unconditionally stable numerical method is based on the operator splitting method (OSM). The nonlinear and linear terms in the gradient flow are solved analytically and using the Fourier spectral method, respectively. The numerical solutions in each step are bounded for any time step size and the overall scheme is temporally second-order accurate. We prove theoretically the unconditional stability and boundedness of the numerical solutions. In addition, several numerical tests are conducted to demonstrate the performance of the proposed method.