Perhaps one way of proceeding is to set this up as a block tridiagonal system of equations, and then solve for all of the unknowns using an efficient algorithm.
I would like to numerically solve the following system of PDEs using a stable computational scheme. The system involves three equations and three unknowns ( , , and ). I've tried to solve this system using an explicit scheme, but I found that the solution could become unstable. An implicit scheme may work, but then how would I deal with the three unknowns in a computationally efficient fashion?
Could someone suggest a stable scheme to numerically solve this system?