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

Plot one trajectory

ClearSystemCache;

varnames := {x1,x2}
a := 1; b := 3; c:=3; d := 1;
odelinear:= {x1'[t] == a x1[t]+b x2[t], x2'[t]== c x1[t]+d x2[t]}; 
x1min := -2 ; x1max := 2; x2min := -2; x2max := 2; 
Tmin :=-1.7;  Tmax  := 1.8;

SolnOde := NDSolve[ {odelinear[[1]], odelinear[[2]], x1[0]== 0.2, x2[0]==0},{x1,x2},{t,Tmin,Tmax}]
ParametricPlot[Evaluate[{x1[t],x2[t]}/.SolnOde],{t,Tmin,Tmax},
PlotRange -> {{x1min,x1max},{x2min,x2max}}, PlotStyle -> Purple, 
Ticks -> {{-2,-1,0,1,2},{-2,-1,0,1,2}}]

Plot several trajectories

odeIC := {{-0.5, 0.5},{-0.5,-0.5},{0.5,-0.5},{0.5,0.5},{ 0.6,0},{-0.6,0},{0,-0.6},{0,0.6}};  
trajectory[t_,x10_,x20_]:=
Block[{x1,x2},
Evaluate[{x1[t],x2[t]}/.NDSolve[{odelinear[[1]],odelinear[[2]],x1[0]==x10,x2[0]==x20},
{x1,x2},{t,Tmin,Tmax}]]]
phaseplot =ParametricPlot[Evaluate[
Table[trajectory[t,iCond[[1]],iCond[[2]]],{iCond,odeIC}] ] , 
{t,Tmin,Tmax},PlotRange->{{x1min,x1max},{x2min,x2max}}, PlotStyle -> Blue,
Ticks->{{-2,-1,0,1,2},{-2,-1,0,1,2}}]

With Direction Field

odevectorfield := {a x1+b x2,c x1+d x2}
range1var := {x1,x1min,x1max}
range2var := {x2,x2min,x2max}
dirgraph=VectorPlot[Evaluate[odevectorfield,range1var,range2var],
VectorScale ->{0.05,Scaled[1.3]},AspectRatio -> 1,Axes -> True,AxesLabel->{range1var[[1]],range2var[[1]]}];
Show[phaseplot,dirgraph]