# Solving a system of PDEs with no-flux boundary conditions (finite difference discretization)

I am interested in solving a system of linear PDEs with the finite difference method and I'm having trouble to solve the no-flux boundary condition correctly. \begin{align} \frac{\partial n}{\partial t} &= D_n\nabla^2n - \nabla\cdot\left(\chi n\nabla c\right) - \nabla\cdot (\rho n\nabla f) \\ \frac{\partial f}{\partial t} &= \beta n-\gamma nf \\ \frac{\partial c}{\partial t} &= -\eta nc \end{align} The domain is $$\Omega=[0,1]\times[0,1]$$ and $$D_n ,\chi, \rho, \beta, \gamma, \eta$$ are constants. Also there is the following boundary condition stated, where $$\xi$$ is the outward unit normal vector: \begin{align} \left[D_n\nabla n - \left(\chi n\nabla c\right) - (\rho n\nabla f)\right]\cdot \xi = 0 \end{align}. Considering $$\xi = (1,0)$$ I received: \begin{align} D_n\partial_xn-(\chi n\partial_xc-\rho n\partial_xf)=0 \Leftrightarrow D_n\partial_xn = \chi n\partial_xc-\rho n\partial_xf \end{align} So now considering on point on the boundary, eg $$(1,n)$$ and discreting with finite difference we get: \begin{align} N_{1,n+1}-N_{1,n-1} = \frac{\chi}{D_n}N_{1,n}(C_{1,n+1}-C_{1,n-1})-\frac{\rho}{D_n}N_{1,n}(F_{1,n+1}-F_{1,n-1}) \\ \Leftrightarrow N_{1,n+1} = N_{1,n-1} + \frac{\chi}{D_n}N_{1,n}(C_{1,n+1}-C_{1,n-1})-\frac{\rho}{D_n}N_{1,n}(F_{1,n+1}-F_{1,n-1}) \end{align} So now we got an equation for our ghost point which we can use for equation $$n$$. The discretization is given as: \begin{align} N_{l,m}^{q+1} = N_{l,m}^{q}P_0 + N_{l+1,m}^{q}P_1+N_{l-1,m}^{q}P_2+N_{l,m+1}^{q}P_3+N_{l,m-1}^{q}P_4 \\ F_{l,m}^{q+1} = F_{l,m}^{q}(1-k\gamma N_{l,m}^{q})+k\beta N_{l,m}^{q} \\ C_{l,m}^{q+1} = C_{l,m}^{q}(1-k\eta N_{l,m}^{q}) \end{align}

if we consider the point $$(1,n)$$ as stated above we'd receive: \begin{align} N_{1,n}^{q+1} = N_{1,n}^{q}P_0 + N_{2,n}^{q}P_1+N_{0,n}^{q}P_2+N_{1,n+1}^{q}P_3+N_{1,n-1}^{q}P_4 \end{align} The index $$(1,n+1)$$ is not in our domain, but we can consider the equation for $$N_{1,n+1}$$ from above and insert. I'm having trouble when considering the $$P_i$$, for example: \begin{align} P_3 = \frac{kD}{h^2}- \frac{k}{4h^2}(\chi(C_{l,m+1}^{q} - C_{l,m-1}^{q})+\rho(F_{l,m+1}^{q}-F_{l,m-1}^{q})) \end{align} My problem is that I don't know how to determine $$C_{1,n+1}$$ and $$F_{1,n+1}$$ to get an equation for $$N_{1,n+1}$$ I can use.

Any suggestions what to try? Please let me know if my question lacks some information, it's my first question here. Any advice is appreciated!