We perform cosmological N-body simulations of the Dvali-Gabadadze-Porrati braneworld model, by solving the full non-linear equations of motion for the scalar degree of freedom in this model, the brane bending mode. While coupling universally to matter, the brane-bending mode has self-interactions that become important as soon as the density field becomes non-linear. These self-interactions lead to a suppression of the field in high-density environments, and restore gravity to General Relativity. The code uses a multi-grid relaxation scheme to solve the non-linear field equation in the quasi-static approximation. We perform simulations of a flat self-accelerating DGP model without cosmological constant. However, the type of non-linear interactions of the brane-bending mode, which are the focus of this study, are generic to a wide class of braneworld cosmologies. The results of the DGP simulations are compared with standard gravity simulations assuming the same expansion history, and with DGP simulations using the linearized equation for the brane bending mode. This allows us to isolate the effects of the non-linear self-couplings of the field which are noticeable already on quasi-linear scales. We present results on the matter power spectrum and the halo mass function, and discuss the behavior of the brane bending mode within cosmological structure formation. We find that, independently of CMB constraints, the self-accelerating DGP model is strongly constrained by current weak lensing and cluster abundance measurements.