# Finite difference scheme for a coupled nonlinear pde

I have a set of two coupled PDEs given by $$I_x^2u+I_xI_yv+I_xI_t=\frac{\partial}{\partial x}\Bigg(\frac{u_x}{\|\nabla u\|}\Bigg)+\frac{\partial}{\partial y}\Bigg(\frac{u_y}{\|\nabla u\|}\Bigg)\\ I_xI_yu+I_y^2v+I_yI_t=\frac{\partial}{\partial x}\Bigg(\frac{v_x}{\|\nabla v\|}\Bigg)+\frac{\partial}{\partial y}\Bigg(\frac{v_y}{\|\nabla v\|}\Bigg)$$ where $I$ is the given input grayscale image and $I_t, I_x$ and $I_y$ are computed by standard methods. From the above equations, I need to compute $u$ and $v$. One way to approach this is to consider the associated evolution equations $$\frac{\partial u}{\partial t}=-(I_x^2u+I_xI_yv+I_xI_t)+\frac{\partial}{\partial x}\Bigg(\frac{u_x}{\|\nabla u\|}\Bigg)+\frac{\partial}{\partial y}\Bigg(\frac{u_y}{\|\nabla u\|}\Bigg)\\\frac{\partial v}{\partial t}=-( I_xI_yu+I_y^2v+I_yI_t)+\frac{\partial}{\partial x}\Bigg(\frac{v_x}{\|\nabla v\|}\Bigg)+\frac{\partial}{\partial y}\Bigg(\frac{v_y}{\|\nabla v\|}\Bigg)$$ I need to solve this iteratively. What would be the correct finite difference approach for this ?

Here is what I have tried in MATLAB:

This is a snippet from the main code. The function mean_curve is here:

My code works only for synthetic images but not for real images. Any help would be appreciated. Thanks.