Phase Portrait of 3d-Systems of ODE
by R Clark Robinson

Plot one trajectory

ClearSystemCache;

varnames := {x1,x2,x3}
odelorenz:= {x1'[t] == sigma (x2[t]-x1[t]), x2'[t]==r x1[t]-x2[t]-x1[t] x3[t],x3'[t]==x1[t] x2[t]-b x3[t]};
sigma := 10; b := 8/3; r := 28;
x1min := -30 ; x1max := 30; x2min := -30; x2max := 30;x3min := 0; x3max := 50; 
Tmin :=0;  Tmax  := 30;
SolnOde := NDSolve[ {odelorenz[[1]],odelorenz[[2]],odelorenz[[3]], x1[0]== 0.1, x2[0]==0.1, x3[0] == 0},{x1,x2,x3},{t,Tmin,Tmax}]
ParametricPlot3D[Evaluate[{x1[t],x2[t],x3[t]}/.SolnOde],{t,Tmin,Tmax},
PlotRange -> {{x1min,x1max},{x2min,x2max},{x3min,x3max}}, PlotStyle -> Blue]

Plot of x1[t] versus t

Plot[Evaluate[x1[t] /. SolnOde], {t,Tmin, Tmax}]

Plot several trajectories

odeIC := {{0.1,0.1,0},{-0.1,-0.1,0}};
trajectory[t_,x10_,x20_,x30_]:=
Block[{x1,x2,x3},
Evaluate[{x1[t],x2[t],x3[t]}/.NDSolve[{odelorenz[[1]],odelorenz[[2]],odelorenz[[3]],x1[0]==x10,x2[0]==x20, x3[0] == x30},
{x1,x2,x3},{t,Tmin,Tmax}]]]
phaseplot =ParametricPlot3D[Evaluate[
Table[trajectory[t,iCond[[1]],iCond[[2]],iCond[[3]]],{iCond,odeIC}] ] , 
{t,Tmin,Tmax},PlotRange->{{x1min,x1max},{x2min,x2max}, {x3min,x3max}}]