The following set of coupled differential + integral equations appear regularly in the literature (e.g Osterbrock, "Astrophysics of Gaseous Nebulae and Active Galactic Nebulae"):
$ n_H(r) \int_{\nu_0}^\infty \frac{L_v}{(4 \pi r^2 h \nu)} (\frac{\nu_0}{\nu})^3 e^{-\tau_\nu}d\nu = Kn_e n_p $
for a known constant K, with ${\nu_0}$ also known.
where we also have the equations
${\tau_\nu} = {\tau_0}(\frac{\nu}{\nu_0})^3 $
$\frac{d\tau_0}{dr} = Mn_H(r) $ (m known constant)
and $L_\nu = L\frac{\nu^3}{e^{k\nu} -1}$ where L and k are also known constants.
$n_e = n_p$
$n_e + n_H(r) = H$, another known constant
In Osterbrock's book he merely states that these equations may be solved computationally and moves on. At this source here http://www.astro.gsu.edu/~crenshaw/4.Ionization.pdf it is claimed "since $\tau_0$ is a known function of $n_H(r)$ and r, the ionisation equilibrium equation can be integrated outwards for given $n_H(r)$ and r to obtain $\frac{n_p}{n_H(r)}$", but I'm not sure I completely understand what is meant by this - we want to solve the integral to find $n_H(r)$ but this itself appears in the integral through tau's dependency on it!
My end goal here is to find $n_H(r)$ as a function of r. My instinct so far is that some substitution should reduce the integral to a purely numerical figure (no $n_H(r)$ dependency) that can be computationally integrated out, but I've had no success with any of the obvious candidates.
Help greatly appreciated!
