h=3; c1=0.2; c2=0.8; c3=3; c4=0.6; 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}];