The set of partial differential equations (1)-(6) is numerically solved by a finite-difference method under adequate initial and boundary conditions.
The numerical schemes used are the same as those described previously [1]}. The methods are based on an explicit-implicit finite difference scheme that is superior in numerical stability. The mesh points \(N_{\rm r}\)  and \(N_{\rm z}\)  in radial and vertical directions are taken to be 410 and 210, respectively, in all simulations.
