Hi, Natchaya,

For instance:

Clear[t];
a = 0.00002963808691235760;
q = 0.00000906440236518953;
b = 0.5;
ts = 120;
system = {S'[t] == (-a*S[t]*Z[t])*(1 + q*Z[t]) + ((b)*Z[t]),
        Z'[t] == (a*S[t]*Z[t])*(1 + q*Z[t]) - ((b)*Z[t])};
initialvalues = {Z[0] == 200, S[0] == 39800};
solution = NDSolve[Join[system, initialvalues], {Z, S}, {t, 0, ts}];
Plot[{Z[t], S[t]} /. solution, {t, 0, ts}, PlotTheme -> "Scientific", 
 PlotRange -> All]
Export["C:\\Users\\Luciano\\Desktop\\Tese\\test.xlsx", 
 Table[Flatten[{t, Z[t], S[t]} /. solution], {t, 0, ts}]]

1. In the first part you must put the values for constants.
2. Then, you write your system (variable "system").
3. In "initialvalues" you put all values to derivatives when t=0.
4. NDSolve will integrate your system. You must inform the range for "t". (For instance: {t, 0, ts})
5. With "Plot" and "export" you can use to generating a figure and file with results.

Best regards.