Here's the deal with Mathematica: in order to cut and paste Mathematica code, I need to first convert it to Raw-Input form which I did and pasted it below but it's hard to read. So cut and paste it back into your Mathematica session, select it, then do a Cell/Convert To/Standard Form so that it's easier to read. I first defined the two equations as eqn1 and eqn2 so make sure I input it correctly. Then defined the parameters and supplied the initial conditions to NDSolve and ran it from t=0 to t=10, obtained the solutions and then did a plot of them which I posted below. Red is q and Blue is z. Also, I don't think you can obtain an explicit solution although Mathematica has DSolve for that. You can try running DSolve[{eqn1,eqn2},{q,z},t]. I tried but stopped it after about one minute as it did not come back with an answer.

Code:

eqn1 = Derivative[1][q][t] ==
(1/(Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + q[t])^2 +
2*Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*(-1 + q[t])*
q[t] - Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s2)*
q[t]^2))*(1 - q[t])*q[t]*(Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*
(-1 + q[t]) + Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*
(-1 + q[t]) + Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*
q[t] - Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s2)*q[t]);
eqn2 = Derivative[1][z][t] == (1/(\[Sigma]^2 + \[Omega]^2)^(3/2))*\[Omega]*
((-Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))])*(-1 + q[t])^2*z[t] -
2*Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*(-1 + q[t])*
q[t]*(a1 + z[t]) + Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*
(-1 + s2)*q[t]^2*(a2 + z[t]));
a1 = 0.5;
a2 = 1;
s1 = 0.0005;
s2 = 0.001;
\[Omega] = 1;
\[Sigma] = 1;
mysol = NDSolve[{eqn1, eqn2, q[0] == 0.001, z[0] == -10}, {q, z},
{t, 0, 10}]
Plot[Evaluate[{q[t], z[t]} /. mysol], {t, 0, 10},
PlotStyle -> {Red, Blue}]