We present an investigation of accumulated numerical error in the simulation of four-wave mixing, which is frequently used to model a variety of nonlinear optical devices. The Dormand-Prince method (commonly used by commercial solving software such as MATLAB) has been found to be susceptible to numerical error, which manifests itself in a 3.8% increase in the total power over 200 m of nonlinear interaction length. This numerical error can lead to qualitatively mistaken physical interpretations of simulation results, which are similar to those found in previously published materials. We use a home-built Adams-Bashforth solver to simulate four-wave mixing, which produces results that do not lead to unphysical results, even for simulation over a large nonlinear interaction length. By comparing the results of these two methods we were able to illustrate the qualitative effects of the accumulated numerical error in the former. The source of this cumulative power error is traced to the solutions provided by the Dormand-Prince method for the self- and cross-phase modulation terms of the coupled mode equations; this error increases for larger nonlinearities or if step size increases. Even when this power accumulates from infinitesimal per-step errors, significant changes occur that could lead to qualitative differences in generated power values and conversion efficiencies. This reveals the potential danger of applying commonly used numerical solvers in simulating nonlinear optical processes.