Phase space for damped pendulum

(* set damping factor *)

b:=0.2;

(* values A,B are the initial conditions for our trajectory *)

Clear[soln]
soln[A_,B_]:=Module[{tmp},
NDSolve[{θ''[t]==-Sin[θ[t]]-b*θ'[t],θ[0]==A,θ'[0]==B},θ,{t,0,25}][[1]]
]

soln[1,0]
θ[1]/.%

(* Plot the fixed points *)

Table[{π*n,0},{n,-4,4,1}]//N
pltFPts=ListPlot[%,PlotMarkers->{Automatic,Medium}]
Cos[%%]

(* Plot trajectories *)

Show[
ParametricPlot[Evaluate[{
{θ[t],θ'[t]}/.soln[2,0],
{θ[t],θ'[t]}/.soln[0,-1],
{θ[t],θ'[t]}/.soln[0,-3],
{θ[t],θ'[t]}/.soln[-4,4],
{θ[t],θ'[t]}/.soln[4,-5],
{θ[t],θ'[t]}/.soln[5,1]
}],{t,0,20},PlotRange->All],pltFPts
]

(* Plot more trajectories *)

Show[
ParametricPlot[Evaluate[{
{θ[t],θ'[t]}/.soln[-1,2.4],
{θ[t],θ'[t]}/.soln[-3.13,0],
{θ[t],θ'[t]}/.soln[3.14,0],
{θ[t],θ'[t]}/.soln[-6,3.45],
{θ[t],θ'[t]}/.soln[6,-3.45],
{θ[t],θ'[t]}/.soln[0,-4],
{θ[t],θ'[t]}/.soln[0,3],
{θ[t],θ'[t]}/.soln[-6,1],
{θ[t],θ'[t]}/.soln[2,1.4],
{θ[t],θ'[t]}/.soln[4,.4]
}],{t,0,20},PlotRange->All]
,pltFPts]