h=3; c1=0.2; c2=0.8; c3=3; c4=0.6; tmax=70; 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=ParametricPlot3D[Evaluate[{m[t],e[t],p[t]}/.loesung], {t,0,tmax},PlotPoints->1000,AxesLabel->{"m","e","p"}, PlotRange->All,AspectRatio->7/10,ImageSize->{400,300}];