The nonlinear propagation and stability of dust ion-waves in plasma is analytically and numerically investigated.By using the the standard reductive perturbation method, the electrostatic potential in dusty pair-ion-electron plasma is modelled by cylindrical Kadomtsev-Petviashvili (CKP) equation. The soliton solutions are obtained using the direct integration for single soliton solution and the Hirota bilinear method to find multisoliton solution of the system. It is noticed that the Hirota method better illustrate the physical reality of dust pair-ion plasma since it generalizes different forms of solutions. From the numerical simulations, it is obseved that, the plasma parameters strongly influence the properties of the soliton solution, namely, the amplitude and the width. The analysis of the stability of the soliton solutions revels that the stable solution copropagates with seven other solutions, eigenmodes of the Legendre equation. These modes contain basic symmetry and axisymmetric configuration consistent with relevant experimental observations in existing experiments.