Phase space

Draw the phase - space for this set of equations, choosing as initial conditions points within the semicircle defined by  x^2+y^2=1 and  y>0 :
 (these are the cosmological equations for a model of coupled quintessence with exponential potential; mu is the exponent of the potential, and beta is the coupling)

mu = λ; beta = Q; (Subscript[x, 1])=x; (Subscript[x, 2])=y; (Subscript[ϱ, r])=0

Everything depends only on x,y: 
omegam=1-x^2-y^2
w=x^2-y^2

xd[x_,y_]:=(-1/2*(3-3*x^2+3*y^2))*x+mu*y^2+beta*(1-x^2-y^2); (* x^2 represents the Kinetic energy of a scalar field *)

yd[x_,y_]:=-mu*x*y-y*(-3-3*x^2+3*y^2)/2;  (* y^2  represents the potential energy of a scalar field  *)

You can also find analytically the fixed points by solving  x’= y’=0, and confirm the graphical result. You can also find the stability of the fixed point by linearizing around the fixed points and finding the eigenvalues of the linear matrix.

tin=-5;tfin=5;
mu=1;beta=0.2;
(* try changing mu, beta *)
sol[xin_,yin_]:=NDSolve[{x'[t]==xd[x[t],y[t]],y'[t]==yd[x[t],y[t]],x[tin]==xin,y[tin]==yin},{x,y},{t,tin,tfin}]

Clear[pp];pp[xin_,yin_]:={xsol=(x/.sol[xin,yin])[[1]];
ysol=(y/.sol[xin,yin])[[1]];ParametricPlot[{xsol[t],ysol[t]},{t,tin,tfin},Frame->True,AspectRatio->1/2,PlotRange->{{-1,1},{0,1}},Axes->False]}[[1]]

Show[Table[pp[xin,.03],{xin,-.99,.99,.05}],Plot[Sqrt[1-x^2],{x,-1,1},PlotStyle->{Thick,Red}]]

Fixed points and stability

Clear[mu,beta];fp=Solve[{xd[x,y]==0,yd[x,y]==0},{x,y}]

list=Chop[fp/.mu->1/.beta->0.2]

?Select

(* select only points in the upper semicircle *)

goodlist=Select[list,#[[2,2]]>=0 && #[[1,2]]<=1&]

(* Find the scalar field EoS for every point: *)

Table[goodlist[[i,1]]^2-goodlist[[i,2]]^2,{i,Length[goodlist]}]

(* linearize around each fixed point *)

linearx=Normal[Series[xd[x+eps*dx,y+eps*dy]/.mu->1/.beta->0.2/.goodlist,{eps,0,1}]]/.eps->1
lineary=Normal[Series[yd[x+eps*dx,y+eps*dy]/.mu->1/.beta->0.2/.goodlist,{eps,0,1}]]/.eps->1

This produces the system 

(dx’,dy’) = M (dx,dy)  

around every critical point

This means that for instance for the first point we have the linearized system

dx' = 3.4*dx + 0*dy
 dy' = 0*dx + 4*dy
 
 so the matrix is {{3.4,0},{0,4}}

(* extract linearization matrix *)
matrix=Table[{{Coefficient[linearx[[i]],dx],Coefficient[linearx[[i]],dy]},{Coefficient[lineary[[i]],dx],Coefficient[lineary[[i]],dy]}},{i,Length[linearx]}]

matrix[[4]]//TableForm

(* find eigenvalues and identify stability *)
eigen=Table[Eigenvalues[matrix[[i]]],{i,Length[linearx]}];
stab=Table[Which[eigen[[i,1]]==0 || eigen[[i,2]]==0,"Degenerate",eigen[[i,1]]>0 && eigen[[i,2]]>0,"Unstable",eigen[[i,1]]*eigen[[i,2]]<0 ,"Saddle",True,"Stable"],{i,Length[linearx]}]

?Labeled

pointlist=Table[{x/.goodlist[[i]],y/.goodlist[[i]]},{i,Length[goodlist]}]

lab=Table[Labeled[pointlist[[i]],stab[[i]],{pointlist[[i]]+{0,-0.04}}],{i,1,4}]

(* Final plot *)

mu=1;beta=0.2;Show[Table[pp[xin,.03],{xin,-.99,.99,.05}],Plot[Sqrt[1-x^2],{x,-1,1},PlotStyle->{Thick,Red}],ListPlot[lab,PlotStyle->{Green,PointSize[.02]}],PlotRange->{{-1.15,1.15},{-.15,1.03}}]