Question: Plotting a three dimensional phase plane

Below is my code for solving a solution of three first order ODEs using the 4th-order Runge-Kutta method.

I have been able to successfully plot the solutions of each of the ODEs (x,y,z) against time t, however I am struggling to produce a plot in the three dimensional phase domain of x,y,z. Could anybody suggest what commands to use as everything I have tried (plot, plot3d, implicitplot3d etc) has produced an error. 

h:= 0.01:
N:= 200:
x:=Vector(N+1):
y:=Vector(N+1):
z:=Vector(N+1):
t:=Vector(N+1):
x[1]:= 0.2:
y[1]:= -0.2:
z[1]:= 0.2:
t[1]:= 0:
for i from 1 to N do
k0:= evalf(h*(-(y[i])*(z[i])+5*(x[i]))): #increment for x
L0:= evalf(h*((x[i])*(z[i])-10*(y[i]))): #increment for y
m0:= evalf(h*((1/3)*(x[i])*(y[i])-3.8*(z[i]))): #increment for z
k1:= evalf(h*(-(y[i]+(L0/2))*(z[i]+(m0/2))+(5*(x[i]+(k0/2))))):
L1:= evalf(h*((x[i]+(k0/2))*(z[i]+(m0/2))-10*(y[i]+(L0/2)))):
m1:= evalf(h*((1/3)*(x[i]+(k0/2))*(y[i]+(L0/2))-3.8*(z[i]+(m0/2)))):
k2:= evalf(h*(-(y[i]+(L1/2))*(z[i]+(m1/2))+5*(x[i]+(k1/2)))):
L2:= evalf(h*((x[i]+(k1/2))*(z[i]+(m1/2))-10*(y[i]+(L1/2)))):
m2:= evalf(h*((1/3)*(x[i]+(k1/2))*(y[i]+(L1/2))-3.8*(z[i]+(m1/2)))):
k3:= evalf(h*(-(y[i]+(L2))*(z[i]+(m2))+5*(x[i]+(k2)))):
L3:= evalf(h*((x[i]+(k2))*(z[i]+(m2))-10*(y[i]+(L2)))):
m3:= evalf(h*((1/3)*(x[i]+(k2))*(y[i]+(L2))-3.8*(z[i]+(m2)))):
x[i+1]:= evalf(x[i]+(1/6)*(k0+(2*k1)+(2*k2)+k3)):
y[i+1]:= evalf(y[i]+(1/6)*(L0+(2*L1)+(2*L2)+L3)):
z[i+1]:= evalf(z[i]+(1/6)*(m0+(2*m1)+(2*m2)+m3)):
t[i+1]:= evalf(t[i]+h):
end do:
with(plots):
plot(t,x,color=blue);
plot(t,y,color=red);
plot(t,z,color=green);