We investigate the dynamics of Gaussian wave packets in multilayer black phosphorus (BP). Time-dependent average position and velocity are calculated analytically and numerically by using a continuum model and a method based on the split-operator technique, respectively. By analyzing the wave-packet trajectories with nonvanishing initial momentum along y direction, we observed transient spatial oscillations due to the effect known as zitterbewegung (ZBW). We also demonstrated, based on the Heisenberg picture and by the calculation of the velocity operators, that the trembling motion along the y direction at small times is unavoidable even for null initial momentum. We verified that the ZBW is directly related to the splitting of the wave packet into two parts moving with opposite velocities, as similar to graphene, and the linear dependence on k_y in out-of-diagonal terms in the Hamiltonian. In addition for phosphorene systems, the two portions of the propagated wave packets have an asymmetric shape for unbalanced ([1,0]^T) and phased different ([1,i]^T)initial pseudospin components, playing an important role in the amplitude, frequency and duration time of the transient oscillations. Electrons in the vicinity of the Fermi energy traveling in N-layer phosphorene also exhibit qualitatively similar trembling motion and transient character as found in the monolayer case, except by the oscillation phase difference and final group velocity achieved after the transient behavior. As a consequence of the anisotropy on the
N-layer BP energy bands, effective masses and group velocities along the x and y directions, the wave packet propagates nonuniformly along the different directions and deforms into an elliptical shape. By comparing our analytical results with those ones obtained by the split-operator technique, we verified a good qualitative and quantitative agreement between them, except for very larger values of wave vector and after long time steps.