PySCo (Python Simulations for Cosmology) runs cosmological N-body simulations across various cosmological models, including Newtonian gravity, MOND gravity (quasi-linear formulation), and parametrized gravity (scale independent). It is a multi-threaded Particle-Mesh code and uses multigrid or fast Fourier transform (FFT) methods. PySCo includes an initial condition generator and power spectrum estimation, and computes the background and growth of density perturbations.