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

Plot one trajectory

ClearSystemCache;

odelinear:= {x1'[t] == a11 x1[t] + a12 x2[t]+a13 x3[t], x2'[t]==a21 x1[t]+a22 x2[t]+a23 x3[t] ,x3'[t]==a31 x1[t] +a32 x2[t]+a33 x3[t]};
a11 := 0; a12 := 3; a13 := 1;
a21 := 4; a22 := 1; a23 := -1;
a31 := 2; a32 := 7; a33 := -5;
x1min := -5 ; x1max := 5; 
x2min := -5; x2max := 5;
x3min :=-5; x3max := 5; 
Tmin :=0;  Tmax  := 10;
SolnOde := NDSolve[ {odelinear[[1]],odelinear[[2]],odelinear[[3]], x1[0]== 1 x2[0]==-1, x3[0] == 5},{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.05,0.05,05},
{-0.05,-0.05,-0.05},
{1.0,-1.0,-5.0},
{-1.0,1.0,5},
{1.05,-0.95,-4.95},
{0.95,-1.05,-5.05},
{-1.01,0.99,4.99},
{-0.99,1.01,5.01}};
trajectory[t_,x10_,x20_,x30_]:=
Block[{x1,x2,x3},
Evaluate[{x1[t],x2[t],x3[t]}/.NDSolve[{odelinear[[1]],odelinear[[2]],odelinear[[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}}]

odeIC2 := {{0.05,0.05,05},
{-0.05,-0.05,-0.05},
{1.05,-0.95,-4.95},
{0.95,-1.05,-5},
{-1.05,0.95,4.95},
{-0.95,1.05,5.00},
{2.05,-1.95,0.05},
{1.95,-2.05,-0.05},
{-2.05,1.95,-0.05},
{-1.95,2.05,0.05}};
phaseplot =ParametricPlot3D[Evaluate[
Table[trajectory[t,iCond[[1]],iCond[[2]],iCond[[3]]],{iCond,odeIC2}] ] , 
{t,Tmin,Tmax},PlotRange->{{x1min,x1max},{x2min,x2max}, {x3min,x3max}}]