Laser Travelling wave rate equations solved numerically

[chitambira]

have a wave equation:

(∂^2 p)/〖∂z〗^2 -1/c^2 (∂^2 p)/〖∂t〗^2

In my case, (fields propagation within a semiconductor laser)
Which can be factored into forward propagating photon density:

(∂p^+)/∂t+C_g (∂p^+)/∂z=C_g Gp^++〖1/2 βR〗_sp

And backward propagation wave (photon density):

〖∂p〗^-/∂t-C_g (∂p^-)/∂t=C_g Gp^-+〖1/2 βR〗_sp

The photons are due to direct modulation of a laser, with the carrier density given by:

∂n/∂t=J/qd-n/τ_nr -〖R_sp-C〗_g G(p^++p^-)

I applied method of lines to discritizze the spatial variable, using backward difference to Forward photon and forward difference to backward photon density to obtain a system of first order ODEs, I then applied first order Finite differences for the time variable to obtain the following equations:

Equation 1
(p^+ (t+1,z)-p^+ (t,z))/Δt+C_g (p^+ (t,z)-p^+ (t-1,z))/Δz=C_g.G〖(n(t,z) ).p〗^++〖1/2 β.R〗_sp (n(t,z))

Equation 2
(p^- (t+1,z)-p^- (t,z))/Δt-C_g (p^- (t,z)-p^- (t,z+1))/Δz=C_g.G〖(n(t,z) ).(p〗^++〖1/2 β.R〗_sp (n(t,z))

Equation 3
(n(t,z)-n(t-1,z))/Δt=J(t,z)/qd-n(t,z)/τ_nr - R_sp (n(t,z) )G〖(n(t,z) )-〖1/2 C〗_g G〖(n(t,z) ) 〖.(p〗〗^+ (t,z)+p〗^- (t,z)+p^+ (t,z-1)+p^- (t,z-1) )

I now want to solve these using MATLAB
Can any body help with any recipe, particularly the best algorithm or code to solve such a scenario varying in both time and space?

I have attached a word document with better readable equations

[Jitse Niesen]

The equations that you wrote down (the ones you get after discretization) are almost the algorithm. They can be translated directly into code. Have you ever programmed a finite difference method before? If not, do you have an example that you can use to get you started?

You don't talk about boundary conditions; you need to think about this before you can implement the method. Also, I notice that you use backward differences for the time variable in the third equation. Is that on purpose or a typo?

[chitambira]

backward difference on third equation is done on purpose, I am trying to use two coupled meshes. I have never programmed a FD before especially in this case where discritization has been done to both the time and spacial variables.
Boundary conditions are that p(0, z) = n(0,z) = 0 (causality)
p^+(t,0) = kp^-(t,0)
p^+(t,L) = kp^-(t,L)

where L, is the total length to the right boundary