Ray tracing through a plane where the refractive index is a function of the distance from the origin Suppose that in a plane parallel to the yz-plane, the index of refraction $n$ is a function of the distance from the origin, $R$, i.e., $n = n(R)$.
We know that (e.g., http://aty.sdsu.edu/explain/atmos_refr/invariant.html), that $n(R)R\sin\theta = \mathrm{constant}$ at every point, where $\theta$ is analogous to the "$z$" angles shown in the diagram:
With the relation above I have the information to calculate $\theta$ at every point, but I'd like to recast it into a differential equation that gives me the path of a ray in terms of $y$ and $z$. Is it possible to find a form for $\frac{dz}{dy}$? I could write, for instance, $\theta = \tan^{-1}\left(\frac{z}{y}\right) - \tan^{-1}\left(\frac{dz}{dy}\right)$. Or would I have to parametrize $y$ and $z$? How would I do this? What is the easiest way of numerically finding the path?
 A: Inspired by the solution given to the brachistochrone problem by the calculus of variations I found a differential equation but in cylindrical plane coordinates $\:(r,\phi)\:$.  Note that the brachistochrone problem is identical to find the light path in a medium with variable refractive index of the form $\:n(y)=\textrm{constant}\cdot y^{-1/2}\:$  where $\:y\:$ a rectilinear coordinate, for example on the vertical line. 
In our case the time to travel from point 1 to point 2 is
\begin{equation}
t_{12}=\int\limits_{1}^{2}\mathrm{d}t=\int\limits_{1}^{2} \dfrac{\mathrm{d}s}{v}
\tag{01}
\end{equation}
where $\:s\:$ the arc-length parameter
\begin{equation}
\mathrm{d}s=\sqrt{\mathrm{d}r^{2}+\left(r\mathrm{d}\phi\right)^{2}}=\sqrt{r^{2}+r'^{\,2}}\:\mathrm{d}\phi,  \quad r'\equiv \dfrac{\mathrm{d}r}{\mathrm{d}\phi}
\tag{02}
\end{equation}
and $\:v\:$ the speed of light, a function of the radius $\:r\:$
\begin{equation}
v(r)=\dfrac{c_{0}}{n(r)}, \quad c_{0}=\textrm{speed of light in empty space}
\tag{03}
\end{equation}
where $\:n(r)\:$ the variable refractive index.
Least time means
\begin{equation}
c_{0}\cdot t_{12}=\int\limits_{1}^{2}L\left(r,r',\phi\right)\mathrm{d}\phi=\text{extremum}
\tag{04}
\end{equation}
where 
\begin{equation}
L\left(r,r',\phi\right) \equiv n(r)\sqrt{r^{2}+r'^{\,2}}
\tag{05}
\end{equation}
is the Lagrangian of the problem.
Note that the Lagrangian does not depend explicitly on the independent variable $\:\phi\:$ so by the Beltrami Identity (1)  the Euler-Lagrange equation is equivalent to this one
\begin{equation}
r'\dfrac{\partial L}{\partial r'}-L = A = \text{constant of the motion}
\tag{06}
\end{equation}
Inserting (05) in (06) we have the differential equation in cylindrical plane coordinates
\begin{equation}
n(r)\dfrac{r^{2}}{\sqrt{r^{2}+r'^{\,2}}}= A = \text{constant}=\textrm{The Refractive Invariant}
\tag{07}
\end{equation}
I believe that comparing this with the Refractive Invariant of the link given we must have (2) 
\begin{equation}
\sin  \mathsf{z} = \dfrac{r}{\sqrt{r^{2}+r'^{\,2}}}= \dfrac{r}{\sqrt{r^{2}+\left(\dfrac{\mathrm{d}r}{\mathrm{d}\phi}\right)^{2}}}
\tag{08}
\end{equation}
The difficulty (or even the possibility) to have  analytical solutions to the differential equation (07) depends upon the function $\:n(r)\:$.
Note that for $\:n(r)=\textrm{constant}\:$, the differential equation (07) becomes
\begin{equation}
\dfrac{r^{2}}{\sqrt{r^{2}+r'^{\,2}}}= A' = \textrm{constant}
\tag{09}
\end{equation}
with a general solution
\begin{equation}
r=\dfrac{A'}{\cos\left(\phi+\phi_{0} \right)}  
\tag{10}
\end{equation}
that is the equation of a straight line as expected (no refraction). 

(1)See here: Beltrami identity
(2)
equation (08) is written as 
\begin{equation}
\sin  \mathsf{z}=\dfrac{1}{\sqrt {1+\left(\dfrac{\mathrm{d}r}{r\mathrm{d}\phi}\right)^{2}}}
\tag{08a}
\end{equation}
that is
\begin{equation}
\cot  \mathsf{z}=\dfrac{\mathrm{d}r}{r\mathrm{d}\phi}
\tag{08b}
\end{equation}
identical to the 3rd equation in Deriving the Differential Equation for Refraction with $\:\phi \equiv \theta, r \equiv R\:$.

We can express the right hand side of equation (08b) in terms of the $\:\left(x,y\right)\:$ cartesian coordinates from equations
\begin{align}
r & =\sqrt{x^{2}+y^{2}}
\tag{09a}\\
\tan \phi & =\dfrac{y}{x}
\tag{09b}
\end{align}
From equations (09)
\begin{align}
\mathrm{d}r & =\dfrac{x\mathrm{d}x+y\mathrm{d}y}{\sqrt{x^{2}+y^{2}}}
\tag{10a}\\
\mathrm{d}\phi & =\dfrac{-y\mathrm{d}x+x\mathrm{d}y}{x^{2}+y^{2}}
\tag{10b}
\end{align}
Inserting (09a),(10a),10b) in (08b) we have 
\begin{equation}
\cot \mathsf{z} =\dfrac{1+\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)\biggl(\dfrac{y}{x}\biggr)}{\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right) - \biggl(\dfrac{y}{x}\biggr)}=\cot \left[\tan^{-1}\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)-\tan^{-1}\biggl(\dfrac{y}{x}\biggr) \right]
\tag{11}
\end{equation}
or
\begin{equation}
\mathsf{z} =\tan^{-1}\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)-\tan^{-1}\biggl(\dfrac{y}{x}\biggr)
\tag{12}
\end{equation}
Replacing $\:\mathsf{z}\rightarrow\theta'\:,\: x\rightarrow y\:,\:y\rightarrow z\:$ we have the OP's equation but with wrong sign
\begin{equation}
\theta' =\tan^{-1}\left(\dfrac{\mathrm{d}z}{\mathrm{d}y}\right)-\tan^{-1}\biggl(\dfrac{z}{y}\biggr)=-\theta
\tag{12}
\end{equation}

(3)  Two Examples :
In the first example, with $\:n(r)=\mathrm{constant}\cdot r^{-2}\:$, the solution of the differential equation (07) is $\: r(\phi)=c_{1}\,\sin(\phi+c_{2})\:$ where $\:c_{1},c_{2}\:$ are constants. This is the parametric equation of a circle. So the light path is a circular arc.

In the second example, with $\:n(r)=\mathrm{constant}\cdot r^{-1}\:$, the solution of the differential equation (07) is $\: r(\phi)=c_{1}\,e^{c_{2}\,\phi} \:$ where $\:c_{1},c_{2}\:$ are positive constants. This is the parametric equation of a spiral. So the light path is a spiral arc.


