Suppose we define the following $\zeta = \ln{\varepsilon}$ and $\xi = \ln{\mu}$, where $\varepsilon$ and $\mu$ are the permittivity and permeability, respectively. In a system with no sources (i.e., $\mathbf{j} = 0$ and $\rho_{c} = 0$), then we know that $\nabla \cdot \mathbf{D} = 0$, where $\mathbf{D} = \varepsilon \ \mathbf{E}$ and $\mathbf{B} = \mu \ \mathbf{H}$. After a little vector calculus we can show that:
$$
\nabla \cdot \mathbf{E} = - \nabla \zeta \cdot \mathbf{E} \tag{0}
$$
Using this and some manipulation of Faraday's law and Ampêre's law, we can show that the general differential equation in terms of electric fields only is given by:
$$
\left( \mu \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \nabla^{2} \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \nabla \zeta + \left( \nabla \zeta \cdot \nabla \right) \mathbf{E} + \nabla \left( \zeta + \xi \right) \times \left( \nabla \times \mathbf{E} \right) \tag{1}
$$
We can get a tiny amount of reprieve from this by assuming that the permeability is that of free space, i.e., $\nabla \xi = 0$. If we further argue that the only direction in which gradients matter is along $\hat{x}$ and that the incident wave vector, $\mathbf{k}$, is parallel to this, then we can further simplify Equation 1 to:
$$
\left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \zeta' \hat{x} + \left( \zeta' \frac{ \partial }{ \partial x } \right) \mathbf{E} + \zeta' \hat{x} \times \left( \nabla \times \mathbf{E} \right) \tag{2}
$$
where $\zeta' = \tfrac{ \partial \zeta }{ \partial x }$.
After some more manipulation, we can break this up into components to show that:
$$
\begin{align}
\text{x : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{x} & = E_{x} \zeta'' + \zeta' \frac{ \partial E_{x} }{ \partial x } \tag{2a} \\
\text{y : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{y} & = 0 \tag{2b} \\
\text{z : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{z} & = 0 \tag{2c}
\end{align}
$$
In the limit where the incident wave is entirely transverse, then $E_{x} = 0$ and the x-component (Equation 2a) is entirely zero.
Next you assume that $\mathbf{E} = \mathbf{E}_{o}\left( x \right) e^{i \omega t}$, where $\omega$ is the frequency of the incident wave. Then there will be incident, reflected, and transmitted contributions to the total field at any given point (well the transmitted is always zero in the first medium, of course). Any incident and transmitted contributions with have $\mathbf{k} \cdot \mathbf{x} > 0$ while reflected waves will satisfy $\mathbf{k} \cdot \mathbf{x} < 0$. You define the ratio of the reflected-to-incident fields (well impedances would be more appropriate) to get the coefficient of reflection.
Simpler Approach
A much simpler approach is to know where to look for the answer to these types of questions. As I mentioned in the comments, there has been a ton of work on this very topic (i.e., spatially dependent index of refraction) done for the ionosphere. If we look at, for instance, Roettger [1980] we find a nice, convenient equation for the reflection coefficient, $R$, as a function of the index of refraction, given by:
$$
R = \int \ dx \ \frac{ 1 }{ 2 \ n\left( x \right) } \frac{ \partial n\left( x \right) }{ \partial x } \ e^{-i \ k \ x} \tag{3}
$$
There is no analytical expression for $R$ for your specific index of refraction. However, numerical integration is not difficult if one knows the values for $d$ and $\epsilon$. Note that if we do a Taylor expansion for small $\epsilon$, then the integrand (not including the exponential) is proportional to cosine, to first order in $\epsilon$ (cosine times sine if we go to second order).
References
- Gossard, E.E. "Refractive index variance and its height distribution in different air masses," Radio Sci. 12(1), pp. 89-105, doi:10.1029/RS012i001p00089, 1977.
- Roettger, J. "Reflection and scattering of VHF radar signals from atmospheric refractivity structures," Radio Sci. 15(2), pp. 259-276, doi:10.1029/RS015i002p00259, 1980.
- Roettger, J. and C.H. Liu "Partial reflection and scattering of VHF radar signals from the clear atmosphere," Geophys. Res. Lett. 5(5), pp. 357-360, doi:10.1029/GL005i005p00357, 1978.
Historically, Cauchy derived his equation assuming that light propagated in an elastic aether. The aether theory of light may be wrong, but Cauchy's derivation was pretty much equivalent to simple models of light propagation in dielectrics. If the resonant absorption frequency of the electrons in the media is $\omega_0$, then the index of refraction $n$ in such a model is given by (e.g. see Feynman Eq. 32.27)
$$n^2=1+\frac{\omega_p^2}{\omega_0^2-\omega^2-i\gamma\omega} $$
where $\omega_p=\sqrt{n_e e^2/m \epsilon_0}$ is the plasma frequency.
As noted in this answer to "Why do we neglect higher order terms in Cauchy's Equation?", in the limit of negligible absorption (i.e. $\gamma=0$), this reduces to
$$n=\sqrt{1+\frac{p^2}{1-x^2}}$$
where $p=\lambda_0/\lambda_p$ and $x=\omega/\omega_0$.
When $\omega<<\omega_0$, this can be expanded in a Taylor series in $x=\omega/\omega_0\approx\lambda_0/\lambda$, which as an even function of $x$ only has even terms in $1/\lambda$:
$$n(\lambda)=\sqrt{p^2 + 1} + \frac{p^2\lambda_0^2}{2 \sqrt{p^2 + 1}} \frac{1}{\lambda^2} + \frac{p^2 (3 p^2 + 4) \lambda_0^4}{8 (p^2 + 1)^{3/2}}\frac{1}{\lambda^4} + O\left(\frac{1}{\lambda^6}\right)$$
This model turns out to be a good enough approximation that empirical fits to $A,B,C$ in
$$ n(\lambda) = A + \frac{B}{\lambda^2} + \frac{C}{\lambda^4} + ...$$
work reasonably well to parameterize $n(\lambda)$ for some common materials for some wavelengths. In particular, it doesn't work badly in the visible for some gases and glasses.
Best Answer
We can rewrite the differential equation as: $$\frac{1}{\lambda}\frac{dn_{eff}}{d\lambda}-\frac{n_{eff}}{\lambda^2}=\frac{d}{d\lambda}\left(\frac{n_{eff}}{\lambda}\right)=-\frac{n_g}{\lambda^2}$$
So$$\frac{n_{eff}}{\lambda}=C-\int\frac{n_g}{\lambda^2}d\lambda$$for some constant $C$, which you can establish using suitable boundary conditions.
To solve this using numerical data, as you have suggested in the comments, between $850\space\text{nm}$ and $855\space\text{nm}$, you would replace the integral with a sum: $$\int\frac{n_g}{\lambda^2}d\lambda\approx\sum\frac{n_g}{\lambda^2}\Delta\lambda$$ where $\Delta\lambda$ is the wavelength spacing between measurements.