My knowledge of finite difference is very basic so this could be very trivial. I've seen how multidimensional finite difference works for say fluid equations, but they are also dealing with a single variable. How would one solve a coupled multidimensional equation, for example of the form
$$a \nabla a = \nabla b,$$
in two dimensions? Meaning we have the system of equations
$$a \, \partial_x a = \partial_x b,$$
$$a \partial_y a = \partial_y b.$$
Using a simplistic finite difference method I can solve each equation independently
$$a(x+\Delta x, y) = \frac{b(x+\Delta x,y) - b(x,y)}{a(x,y)} + a(x,y),$$
and
$$a(x, y+\Delta y) = \frac{b(x, y+\Delta y) - b(x,y)}{a(x,y)} + a(x,y).$$
However if I want to get to $a(x+\Delta x, y+\Delta y)$, there are two options either first moving along $x$ and then $y$ or vice versa. In this case I am not guaranteed to get the same result, meaning the answer is path dependent which is not desirable. What is the common practice to deal with this?
Edit: I'll compute the values over a square grid with initial conditions $a(0,0) = a_0$ and $b(x,y)=x-y$.
My purpose is to solve this equation so that I may initialize the grid in $a$. This will be implemented in a program I am writing. The equation given is just a no-frills re-expression of what I'm attempting to do. My true equation is solving an isentropic hydrostatic atmosphere with a complicated potential, $\nabla \Phi$, for the density structure, $\rho$. The equation reduces to $$\gamma \, K \, \rho^{\gamma-2} \, \vec{\nabla} \rho = - \vec{\nabla} \Phi$$, where $K$ and $\gamma$ are some constants.