The instabilities observed in direct numerical simulations of the Gross–Neveu equation under linear and harmonic potentials are studied. The Lakoba algorithm, based on the method of characteristics, is performed to numerically obtain the two spinor components. We identify non-conservation of energy and charge in simulations with instabilities, and we find that all studied solitons are numerically stable, except the low-frequency solitons oscillating in the harmonic potential over long periods of time. These instabilities, as in the case of the Gross–Neveu equation without potential, can be removed by imposing absorbing boundary conditions. The dynamics of the soliton is in perfect agreement with the prediction obtained using an Ansatz with only two collective coordinates, namely, the position and momentum of the center of mass. We employ the temporal variation of both field energy and momentum to determine the evolution equations satisfied by the collective coordinates. By applying the same methodology, we also demonstrate the spurious character of the reported instabilities in the Alexeeva–Barashenkov–Saxena model under external potentials.