This paper presents a numerically stable, time-accurate algorithm for simulating the gradient flow associated with the Modica–Mortola functional with a uniformly spaced multi-well potential. The scheme uses operator splitting; the nonlinear component is updated analytically, while the linear part is advanced by a Fourier spectral discretization. The method is unconditionally stable, preserves pointwise boundedness independently of the time step size, and attains spectral accuracy in space and first-order accuracy in time. We provide a theoretical analysis establishing unconditional stability and boundedness, and present comprehensive numerical experiments that demonstrate the accuracy and robustness of the proposed approach.