I am trying to reproduce the analysis in the paper, "Polariton condensation threshold investigation through the numerical resolution of the generalized Gross-Pitaevskii equation". I have a doubt about how the time dependence enters in the RHS of Eqn 7 and 8 when one is using the RK4 algorithm to solve the coupled equations. I would greatly appreciate it if someone could provide a pseudocode to solve this system of equations using the finite-difference and RK4 methods.