mu[s_] := M*s/(a + s)

M = 
      10

a = 0.5

mu[s]

dil[t_] := 
      Sin[t]/(2 + Cos[t])  +  
        M*(2 - Cos[t])/(4*a  +  2  - 
                Cos[t])

dil[t]

graf1 = Plot[dil[t], {t, 0, 30}]

Display["dilution-new", graf1, "EPS"]

xref[t_] := 0.5 + 0.25*Cos[t]

xref[t]

nutin[t_] := 1

nutin[t]

Plot[nutin[t], {t, 0, 50}]

result = 
      NDSolve[{
z'[t] ==  
            dil[t]*(nutin[t] - z[t]), x'[
              t]  ==  (mu[z[t] - x[t]] -  dil[t])*x[t], 
          z[0] == 3, x[0] == 2}, {z, x}, {t, 0, 100}]

graf3 = Plot[x[t] /. result, {t, 0, 30}];

Display["x-solution-new", graf3, "EPS"]

graf4 = Plot[xref[t], {t, 0,  30}]

Display["x-ref-new", graf4, "EPS"]

graf5 = Plot[z[t] /. result, {t, 0, 10}];

Display["z-solution-new", graf5, "EPS"]