I am trying to solve the following differential equation describing the steady-state carrier density and photon density in a laser diode: $$AN(z)+BN(z)^2+CN(z)^3=\frac{\eta_iI}{qV}-v_gg(N)(S^+(z)+S^-(z))$$ $$\frac{dS^+(z)}{dz}=(\Gamma g(N)-\alpha_i)S^+(z)+\frac{\Gamma R_{sp}'}{v_g}$$ $$-\frac{dS^-(z)}{dz}=(\Gamma g(N)-\alpha_i)S^-(z)+\frac{\Gamma R_{sp}'}{v_g}$$ The reflections at the facets of the laser define the following boundary conditions: $$S^+(0)=R_1S^-(0)$$ $$S^-(L)=R_2S^+(L)$$ I am trying to use the finite-difference method, where I divide the laser into $M+1$ points, so I make the following approximation: $$\frac{dS}{dz}\approx\frac{S(z+\Delta z)-S(z-\Delta z)}{2\Delta z}\equiv\frac{S_{i+1}-S_{i-1}}{2\Delta z}$$ This gives me a system of $3(M+1)$ equations: $$AN_i+BN_i^2+CN_i^3=\frac{\eta_iI}{qV}-v_gg(N_i)(S_i^++S_i^-),\quad(i=1,...,M+1)$$ $$\frac{S_{i+1}^+-S_{i-1}^+}{2\Delta z}=(\Gamma g(N_i)-\alpha_i)S_i^++\frac{\Gamma R_{sp}'}{v_g},\quad(i=2,...,M)$$ $$-\frac{S_{i+1}^--S_{i-1}^-}{2\Delta z}=(\Gamma g(N_i)-\alpha_i)S_i^-+\frac{\Gamma R_{sp}'}{v_g},\quad(i=2,...,M)$$ $$S^+_{1}=R_1S^-_{1}$$ $$S^-_{M+1}=R_2S^+_{M+1}$$ This is where I run into a problem. I have two more equations I need to define: $$\frac{S_{M+2}^+-S_{M}^+}{2\Delta z}=(\Gamma g(N_{M+1})-\alpha_i)S_{M+1}^++\frac{\Gamma R_{sp}'}{v_g}$$ $$\frac{S_{2}^--S_{-1}^-}{2\Delta z}=(\Gamma g(N_{1})-\alpha_i)S_{1}^-+\frac{\Gamma R_{sp}'}{v_g}$$ However, this leads to two variables that are "ill-defined", namely $S_{-1}^-$ and $S_{M+2}^+$.
I tried to solve this problem by using forward and backward differences at the edges, instead of central differences, to get rid of these out-of-the-boundary variables. However, this reduces the accuracy of the solution significantly. In fact, for $M=100$, I will then get the following result for $N$, $S^+$, $S^-$:
I got this result by using Matlab's symbolic numerical solver vpasolve. Changing the initial values did not change the result, so I suspect that the result is stable. It also satisfies the boundary conditions nicely. However, the values of $N$ at the edges are not physical, and especially $N$ (in the middle) and $S^+$ have pronounced oscillation, which should not exist in the real laser. I suspect this is due to the fact that I use the less-accurate forward and backward differences at the edges, which deteriorates the accuracy of the whole system. I could of course increase the number of points, but the calculation is heavy as it is.
I would like to use central differences at the edges as well, but I am not sure how it works, because of the two additional variables $S_{-1}^-$ and $S_{M+2}^+$, as I am not sure how to handle them. This would also mean that I would have $3(M+1)+2$ variables, but only $3(M+1)$ equations, so would I not then get multiple possible solutions?
Any hints are appreciated.
Thank you in advance.