A complete treatment of your question (more than you ever wanted) is given at http://homepage.tudelft.nl/q1d90/FBweb/diss.pdf, especially section 2.1.1 "Differential equation of light rays in inhomogeneous media". Trying to extract the most useful expression from that dissertation, I believe that the equation you are looking for is:

$$\nabla \Phi a = \frac{2\pi}{\lambda}n(R)$$

where
$\Phi$ = phase
$R$ = position vector
$a$ = unit vector pointing along ray

With a bit of manipulation, that turns into

$$\frac{d}{ds}\left(n\frac{dR}{ds}\right) = \nabla n$$ (equation 2.1.8 in the above reference).

The factor $ds$ can be a bit tricky since it is pointing along the ray - if you want things in X,Y coordinates then you need to worry about the length of $ds$ when it is no longer at a small angle to the X axis - it becomes $\sqrt{dx^2+dy^2}$

A complete treatment of your question (more than you ever wanted) is given here, especially section 2.1.1 "Differential equation of light rays in inhomogeneous media". Trying to extract the most useful expression from that dissertation, I believe that the equation you are looking for is:

$$\nabla \Phi a = \frac{2\pi}{\lambda}n(R)$$

where
$\Phi$ = phase
$R$ = position vector
$a$ = unit vector pointing along ray

With a bit of manipulation, that turns into

$$\frac{d}{ds}\left(n\frac{dR}{ds}\right) = \nabla n$$ (equation 2.1.8 in the above reference).

The factor $ds$ can be a bit tricky since it is pointing along the ray - if you want things in X,Y coordinates then you need to worry about the length of $ds$ when it is no longer at a small angle to the X axis - it becomes $\sqrt{dx^2+dy^2}$