h=3; c1=0.2; c2=0.8; c3=3; c4=0.1; tmax=35;
m0=1; e0=1; p0=1;
modell={m'[t]==1/(c1+p[t]^h) - m[t],
        e'[t]== m[t] - c2 e[t],
        p'[t]== e[t] - c3 p[t]/(c4+p[t]),
        m[0]==m0, e[0]==e0, p[0]==p0};
loesung=NDSolve[modell,{m,e,p},{t,0,tmax},MaxSteps->2000];
bild=Plot[Evaluate[{m[t],p[t]}/.loesung],{t,0,tmax},PlotRange->All,
AxesLabel->{"t","m,p"},
PlotStyle->{{},{Dashing[{0.01}]}},ImageSize->{400,300}];