Cálculo numérico de la generación del segundo armónico en medios no lineales

Recientemente he estado trabajando en el cálculo numérico de la generación de segundos armónicos en medios no lineales. Me refiero al libro Laser Beam Propagation in Nonlinear Optical Media escrito por Shekhar Guha y Leonel P. Gonzalez, CRC Press 2013.

En el capítulo 7 del libro, propusieron un algoritmo de cálculo para resolver las siguientes ecuaciones.

$$\begin{align} \frac{\partial A_p}{\partial z} &= \frac{i}{2k_p} \nabla^2_{T}A_p + \frac{i2d_\text{eff} \omega_p}{cn_p}A^*_p A_s e^{i(K_s - 2K_p)z} -\frac{\alpha_p}{2}A_p \\ \frac{\partial A_s}{\partial z} &= \frac{i}{2k_s} \nabla^2_{T}A_s + \frac{i2d_\text{eff} \omega_p}{cn_s}A_p A_p e^{-i(K_s - 2K_p)z} -\frac{\alpha_s}{2}A_p \end{align}$$ Simplemente simplifique las ecuaciones anteriores usando operadores lineales y no lineales $$\begin{align} \frac{\partial A_p}{\partial z} &= \hat{P_p}A_p + \hat{NL_p} \\ \frac{\partial A_s}{\partial z} &= \hat{P_s}A_s + \hat{NL_s} \end{align}$$ donde estan los operadores $$\begin{align} \hat{P_p} &= \frac{i}{2k_p} \nabla^2_{T} \\ \hat{NL_p} &= \frac{i2d_\text{eff} \omega_p}{cn_p}A^*_p A_s e^{i(K_s - 2K_p)z} -\frac{\alpha_p}{2}A_p \\ \hat{P_s} &= \frac{i}{2k_s} \nabla^2_{T} \\ \hat{NL_s} &= \frac{i2d_\text{eff} \omega_p}{cn_s}A_p A_p e^{-i(K_s - 2K_p)z} -\frac{\alpha_s}{2}A_p \end{align}$$ Entonces le corresponde al método de paso dividido resolver las ecuaciones.

Divida la propagación en N rebanadas, cada una de longitud$\Delta z$

1. Determinar el campo de bombeo incidente$A_p (x, y, z = 0, t)$
2. Colocar$A_s(x,y,z=0,t) = 0$
3. Colocar$\hat{NL} = 0$y propagar los campos usando$\hat{P}$por una distancia$\Delta z$. De$A_s(x,y,z=j\Delta z,t)$Llegar$A_s(x,y,z=(j+1)\Delta z,t)$, podemos usar la transformada casi rápida de Hankel para obtener los resultados.
4. Colocar$\hat{P} = 0$y propagar los campos usando$\hat{NL}$por una distancia$\Delta z$utilizando técnicas de diferencias finitas (discutidas a continuación) para determinar$A_P$y$A_S$.
5. El$A_P$y$A_S$campos se convierten en las entradas para el siguiente$\Delta z$rebanada.
6. Repita 3, 4 y 5 hasta que los campos se propaguen hasta el final del cristal.
7. Calcule los campos fuera del cristal usando el coeficiente de transmisión apropiado.

Ahora viene la pregunta:

¿Debo usar el$A_s(x,y,(j+1)\Delta z,t)$y$A_p(x,y,(j+1)\Delta z,t)$que calculé en el paso 3, o en el paso 4?

El detalle del paso 4 descrito en el libro se enumera aquí: Configuración$\hat{P} = 0$, resuelve la parte no lineal de las ecuaciones.

$$\begin{align} \frac{\partial A_p}{\partial z} &= \frac{i2d_\text{eff} \omega_p}{cn_p}A^*_p A_s e^{i(K_s - 2K_p)z} -\frac{\alpha_p}{2}A_p \\ \frac{\partial A_s}{\partial z} &= \frac{i2d_\text{eff} \omega_p}{cn_s}A_p A_p e^{-i(K_s - 2K_p)z} -\frac{\alpha_s}{2}A_p \end{align}$$

Estas ecuaciones se resuelven utilizando el método implícito y directo de diferencias finitas a partir de los valores conocidos expresados ​​en forma de diferencias finitas como

$$\begin{align} \frac{A_s^{(j+1)\Delta z} - A_s^{(j)\Delta z}}{\Delta z} &= C_s A^{(j+1)\Delta z}_p A_p^{(j+1)\Delta z} -\frac{\alpha_s}{2}A^{(j+1)\Delta z}_s \\ \frac{A_p^{(j+1)\Delta z} - A_p^{(j)\Delta z}}{\Delta z} &= C_s A^{*,j\Delta z}_p A_s^{(j+1)\Delta z} -\frac{\alpha_p}{2}A^{(j+1)\Delta z}_p \end{align}$$

dónde

$$\begin{align} C_p &= \frac{i2d_\text{eff} \omega_p}{cn_p} e^{i(K_s - 2K_p)(j \Delta z)} \\ C_s &= \frac{i2d_\text{eff} \omega_p}{cn_s} e^{-i(K_s - 2K_p)(j \Delta z)} \end{align}$$

Me confundió que ya obtuve el resultado calculado de la$A_s(x,y,(j+1)\Delta z,t)$y$A_p(x,y,(j+1)\Delta z,t)$hasta el paso 3. Y tengo que calcular el$A_s(x,y,(j+1)\Delta z,t)$y$A_p(x,y,(j+1)\Delta z,t)$de nuevo en el paso 4 utilizando el resultado de la$A_s(x,y,(j+1)\Delta z,t)$y$A_p(x,y,(j+1)\Delta z,t)$que ya se calculó en el paso 3.

¿Alguien puede ayudarme con este problema? ¿Cuál es la relación entre el paso 3 y el paso 4?