We present an analytic theory of inhomogeneous superconducting pairing in strongly disordered materials, which are moderately close to Superconducting-Insulator Transition. Within our model, single-electron eigenstates are assumed to be Anderson-localized, with a large localization volume. Superconductivity then develops due to coherent delocalization of originally localized preformed Cooper pairs. The key assumption of the theory is that each such pair is coupled to a large number Z >> 1 of similar neighboring pairs. We derived integral equations for the probability distribution P\(\Delta) of local superconducting order parameter \Delta(r) and solved them analytically in the limit of small dimensionless Cooper coupling constant \lambda << 1. The shape of the order-parameter distribution is found to depend crucially upon the effective number of "nearest neighbors", expressed as the number of neighbors with single-particle energies being within the strip of width 2 \Delta near the given state. The solution we provide is valid for both large and small effective numbers of neighbors. One of our key findings is the discovery of a broad range of parameters where the distribution function P(\Delta) is non-Gaussian but also free of "fat tails" and other features of criticality. The analytic results are supplemented by numerical data, and good agreement between them is observed.