Stochastic epidemic models on networks are inherently high-dimensional and the resulting exact models are intractable numerically even for modest network sizes. Mean-field models provide an alternative but can only capture average quantities, thus offering little or no information about variability in the outcome of the exact process. In this paper we conjecture and numerically prove that it is possible to construct PDE-limits of the exact stochastic SIS epidemics on regular and Erdős-Rényi networks. To do this we first approximate the exact stochastic process at population level by a Birth-and-Death process (BD) (with a state space of $O(N)$ rather than $O(2^N)$) whose coefficients are determined numerically from Gillespie simulations of the exact epidemic on explicit networks. We numerically demonstrate that the coefficients of the resulting BD process are density-dependent, a crucial condition for the existence of a PDE limit. Extensive numerical tests for Regular and Erdős-Rényi networks show excellent agreement between the outcome of simulations and the numerical solution of the Fokker-Planck equations. Apart from a significant reduction in dimensionality, the PDE also provides the means to derive the epidemic outbreak threshold linking network and disease dynamics parameters, albeit in an implicit way. Perhaps more importantly, it enables the formulation and numerical evaluation of likelihoods for epidemic and network inference as illustrated in a worked out example.