This paper presents a further numerical study of the interaction dynamics for solitary waves in a nonlinear Dirac model with scalar self-interaction, the Soler model, by using a fourth order accurate Runge-Kutta discontinuous Galerkin method. The phase plane method is employed for the first time to analyze the interaction of Dirac solitary waves and reveals that the relative phase of those waves may vary with the interaction. In general, the interaction of Dirac solitary waves depends on the initial phase shift. If two equal solitary waves are in-phase or out-of-phase initially, so are they during the interaction; if the initial phase shift is far away from 0 and π, the relative phase begins to periodically evolve after a finite time. In the interaction of out-of-phase Dirac solitary waves, we can observe: (a) full repulsion in binary and ternary collisions, depending on the distance between initial waves; (b) repulsing first, attracting afterwards, and then collapse in binary and ternary collisions of initially resting two-humped waves; (c) one-overlap interaction and two-overlap interaction in ternary collisions of initially resting waves.