If I recall correctly (from about a year ago), there is an interaction between maxmesh and abserr. Depending on the difficulty of the problem, you can go only so far in getting accuracy.

Needs a lot of work, but something that generates the legend separately, such as with(plots): p1:=plot3d(phi*cos(phi)^3,theta=0..2*Pi,phi=0..Pi/2, coords=spherical,scaling=constrained,color=phi,style=patchnogrid,axes=normal): p2:=plot3d(phi,phi=0..Pi/2,y=0..1,color=phi,axes=boxed,style=patchnogrid,orientation=[180,90]): display(p1);display(p2); This doesn't put them together very well because one is in spherical, the other in cartesian coordinates. With quite a bit of work you could get them together on the same plot and co-ordinates. The colors may not correspond here because of the automatic scaling. You may need to look at the colorfunc information and write your own procedure. Maple is great for looking at stuff, but not so good at prettying it up.

Needs a lot of work, but something that generates the legend separately, such as with(plots): p1:=plot3d(phi*cos(phi)^3,theta=0..2*Pi,phi=0..Pi/2, coords=spherical,scaling=constrained,color=phi,style=patchnogrid,axes=normal): p2:=plot3d(phi,phi=0..Pi/2,y=0..1,color=phi,axes=boxed,style=patchnogrid,orientation=[180,90]): display(p1);display(p2); This doesn't put them together very well because one is in spherical, the other in cartesian coordinates. With quite a bit of work you could get them together on the same plot and co-ordinates. The colors may not correspond here because of the automatic scaling. You may need to look at the colorfunc information and write your own procedure. Maple is great for looking at stuff, but not so good at prettying it up.

Just add 'maxmesh'=500 as an option in dsolve. Documented under ?dsolve,bvp