Collective effects with strong coupling between the longitudinal and transverse dynamics are of fundamental importance for applications of high-intensity bunched beams. The self-consistent Vlasov–Maxwell equations are applied to high-intensity bunched beams, and a generalized δ f particle simulation algorithm is developed for bunched beams with or without energy anisotropy. Numerically, the distribution function is spit into a reference distribution and a perturbed part. The perturbed distribution function is represented as a weighted summation over discrete particles, where the particle orbits are advanced by the equations of motion in the focusing field and self-generated fields, and the particle weights are advanced by an equation equivalent to the nonlinear Vlasov equation. The nonlinear δ f method exhibits minimal noise and accuracy problems in comparison with standard particle-in-cell simulations. Systematic studies are carried out for the particle dynamics under conditions corresponding to strong 3D nonlinear space-charge force. The simulations showed that finite bunch-length effects on the collective excitations become insignificant when the aspect ratio ( z b / r b ) is larger than 10 for a moderately intense beam with normalized intensity s b = ω p b 2 / 2 ω β 2 = 0.27 . For bunched beams with energy anisotropy ( T ∥ / T ⊥ < 1 ) , a reference state has been constructed and a dynamic equilibrium is established in the simulations. Collective excitations relative to the dynamic equilibrium have also been successfully simulated by the generalized δ f algorithm.