I am plotting a solution for a known chaotic system of three first order ODEs, along with their isosurfaces, which are obtained by setting the derivative equal to zero for each of the three equations. The change in the solution in the x, y and z direction should be zero through the surfaces Isurf1, Isurf2 and Isurf3 respectively. When plotted, this displays (usually) incorrectly: DE1:=diff(x(t),t)=-z(t); DE2:=diff(y(t),t)=x(t)-y(t); DE3:=diff(z(t),t)=alpha*x(t)+y(t)^2+beta*z(t); DE1a:=subs(alpha=3.1,beta=0.5,DE1); DE2a:=subs(alpha=3.1,beta=0.5,DE2); DE3a:=subs(alpha=3.1,beta=0.5,DE3); SPROTT:=DEplot3d({DE1a,DE2a,DE3a},{x(t),y(t),z(t)}, t=0..300,stepsize=.05,[[x(0)=.1,y(0)=.1,z(0)=.1]],linecolor=black,thickness=1): Isurf1:=implicitplot3d({z=0}, x=-9..9, y=-9..9, z=-9..9, axes=boxed,color=blue): Isurf2:=implicitplot3d({x-y}, x=-9..9, y=-9..9, z=-9..9,axes=boxed,color=yellow): Isurf3:=implicitplot3d({z=-6.2*x-2*y^2}, x=-9..9, y=-9..9, z=-9..9, axes=boxed,color=red): display(Isurf1,Isurf2,Isurf3,SPROTT); However, when the variables y and z are interchanged in the definitions of the ODEs, DE1:=diff(x(t),t)=-y(t); DE2:=diff(z(t),t)=x(t)-z(t); DE3:=diff(y(t),t)=alpha*x(t)+z(t)^2+beta*y(t); the plot displays correctly (usually). What is especially troubling is that sometimes it plots correctly and sometimes it doesn't. I have done this literally dozens of times and can find no pattern or clue as to when I will need to interchange the variables y and z to get this to plot properly.