I am faced with the following system of coupled nonlinear partial differential equations $$ \begin{array}{lcccl} \varphi_{tt} &-& a_1\varphi_{xx} &+& a_2\varphi_t &=& a_3\sin{\varphi} + a_4\sin{\left(\frac{\varphi}{2}\right)} + a_5u_x\sin{\varphi} + a_6v_x\cos{\varphi} + a_7\sin{(\omega t)},\\ u_{tt} &-& b_1u_{xx} &+& b_2u_t &=& b_3(\cos{\varphi})_x, \\ v_{tt} &-& c_1v_{xx} &+& c_2v_t &=& c_3(\sin{\varphi})_x, \end{array} $$ where $a_k,b_k,c_k\in\mathbb{R},~a_1,b_1,c_1>0.$
What ideas might here be for numerically solving this system?
UPD1. Regarding to Chris' question.
Parameter $\omega$ is also constant, same as $a_k,b_k,c_k.$ Indeed, $\left(\cos{\varphi}\right)_x=-\varphi_x\sin{\varphi}.$
I thought to numerically investigate the soliton propagation, and the initial and boundary conditions may follow from here.
I haven't written down the difference approximation yet because, perhaps, it's not the central differences that need to be applied, but upwind-type schemes that I don't know yet.