We develop a theory of the local density of states (LDOS) of disordered superconductors, employing the non-linear sigma-model formalism and the renormalization-group framework. The theory takes into account the interplay of disorder and interaction couplings in all channels, treating the systems with short-range and Coulomb interactions on equal footing. We explore 2D systems that would be Anderson insulators in the absence of interaction and 2D or 3D systems that undergo Anderson transition in the absence of interaction. We evaluate both the average tunneling density of states and its mesoscopic fluctuations which are related to the LDOS multifractality in normal disordered systems. The obtained average LDOS shows a pronounced depletion around the Fermi energy, both in the metallic phase (i.e., above the superconducting critical temperature Tc) and in the insulating phase near the superconductor-insulator transition (SIT). The fluctuations of the LDOS are found to be particularly strong for the case of short-range interactions -especially, in the regime when Tc is enhanced by Anderson localization. On the other hand, the long-range Coulomb repulsion reduces the mesoscopic LDOS fluctuations. However, also in a model with Coulomb interaction, the fluctuations become strong when the systems approaches the SIT.