I need to solve the following system of equations for subject to the given initial and boundary conditions. Please give me some idea which numerical scheme I should use for better accuracy and how should I proceed for the numerical solution of this equation. The coupled boundary conditions are challenging for me. Please help.

By using Crank Nicholson method and Thomas algorithm I can solve the following, but the above one gives trouble for me.