The piecewise-linear appearance of the plots can be removed by using **odeplot **with option **refine**. This will greatly reduce the number of points needed for smoothness, and it will figure out the number of variable-stepsize points on its own. It's an implementation of something akin to the **adaptive** plotting of the regular **plot **command.

All of your plots *can be done* and *should be done* with **odeplot**. The **refine **option must be used in conjunction with the **range** option in the **dsolve** command. It cannot be used with the **numpoints** option.

If you're going to do multiple plots over different ranges, then, of course, the **range **you pass to **dsolve** should be a superset of them all. I changed the order of your plots below to avoid an ugly-looking warning. The warning happens when a plot with a small range is followed by a plot with a larger range. The warning can be completely ignored if the larger range fits in the **range** passed to **dsolve**. I suspect that the warning message was put in the code before the **refine** option existed.

Several other tips:

- I relaxed the error tolerances to
**1e-5** from your **1e-8**. This is sufficient for accurate plotting, and it reduced the number of points in the largest plots by a factor of 4.
- Use the
**parameters** option to **dsolve**. I've shown you this before.
- "Length of output exceeds ..." is not an error message or even a warning. It's an notification that Maple has sucessfully completed a calculation, but won't display it results because the amount of displayed data could seriously degrade the performance of the GUI display. The results of the calculation are still fully accessible.

MaplePrimes will not let me display this worksheet inline, I guess due to the large plots. You will see in the worksheet that they are all smooth.

Download OdeplotRefine.mw

Here is the complete code of the worksheet:

**restart:
sys:= {
diff(x(t),t) = x(t)*(r - b*x(t) - y(t)*(c + beta/(a+x(t)))),
diff(y(t),t) = y(t)*(beta*x(t)/(a+x(t)) - mu)
}:
ics:= {(x,y)(0)=~ (0.2, 0.05)}:
Param:= [r,b,c,a,mu,beta]:
ParamV:= [1, 1, 0.01, 0.36, 0.4, 0.75]:
sol:= dsolve(
sys union ics, {x,y}(t), parameters= Param, numeric, maxfun= 0,
range= 0..20000, (abserr, relerr)=~ 1e-5
):
sol(parameters= ParamV);
plots:-odeplot(sol, [t, (x,y)(t)], t= 0..20000, refine= 1);
#Number of computed points in plot:
Npts:= P-> [rtable_dims]~([indets(P, rtable)[]]):
Npts(%%);
plots:-odeplot(
sol, `[]`~(t, [x,y](t)), t= 0..400, axes= boxed, gridlines,
color= ["Blue", "Red"], title= "Trajectories\n", legend= [x,y](t)
);
Npts(%);
eq_points_sym:= solve(eval(sys, diff= 0), {x,y}(t));
eq_points:= eval([eq_points_sym], Param=~ ParamV);
plots:-odeplot(
sol, [x,y](t), t= 0..20000, refine= 1, color= red, thickness= 2,
axes= boxed, gridlines, title= "Phase-space Diagram\n"
);
Npts(%);
#Task 2
ParamV:= [1, 1, 0.01, 0.36, 0.4, 1]:
sol(parameters= ParamV);
plots:-odeplot(sol, [t, (x,y)(t)], t= 0..20000, refine= 1);
Npts(%);
plots:-odeplot(
sol, `[]`~(t, [x,y](t)), t= 0..400, axes= boxed, gridlines, refine= 1,
color= ["Blue", "Red"], title= "Trajectories\n", legend= [x,y](t)
);
Npts(%);
plots:-odeplot(
sol, [x,y](t), t= 0..20000, refine= 1, color= red, thickness= 0.9,
axes= boxed, gridlines, title= "Phase-space Diagram\n"
);
Npts(%);
eq_points:= eval([eq_points_sym], Param=~ ParamV);
**

And here is my new version of the piecewisey plot that you showed in the Question. This is done with only 8733 points.