I am trying to implement an angionesis model described by Anderson and Chaplin in 1998. The model is based on a set of PDEs defined on an unit square with the following no-flux boundary condition defined on the boundaries of the square:
$$\mathbf{n}\cdot \left(-D\nabla h + h\left(\frac{\chi}{1+\alpha c}\nabla c+\rho\nabla f\right)\right)=0$$ where $\mathbf{n}$ is an outward unit normal vector, $h$, $c$ and $f$ are the unknown functions and $D$, $\chi$, $\alpha$ and $\rho$ are constant parameters.
The problem is that when trying to impose this boundary condition in a finite difference method scheme for $c$ or $f$ functions the function $h$ will be dividing. For example, for the right boundary, $\mathbf{n} = (1,0)$, the equation will be
$$\frac{\partial f}{\partial x} = \frac{1}{h\rho}\left(D\frac{\partial h}{\partial x}-h\left(\frac{\chi}{1+\alpha c}\frac{\partial c}{\partial x}\right)\right),$$ and after discretization and using ghost points the boundary condition can be inposed in the problem. However, the unknown function $h$ is usually equal to zero in most parts of the domain so problems will arise because it is dividing.
How do you usually deal with this type of problem? Any particular approach?
Thanks!