(*Generate blue noise to sample the plane*)
range = 5;
blue = {RandomReal[{-range, range}, {2}]};
Do[
  n = Length[blue];
  candidates = RandomReal[{-range, range}, {n + 1, 2}];
  bestcandidatepos = 
   Position[
     Table[Min[Norm[candidates[[j]] - #] & /@ blue], {j, 1, n}], Max[Table[Min[Norm[candidates[[j]] - #] & /@ blue], {j, 1, n}]] ][[1, 1]];
  AppendTo[blue, candidates[[bestcandidatepos]]];
  , 10^2];
(*definitions*)
eqs[matrix_] := ({q'[t], p'[t]} == matrix . {q[t], p[t]});
initialpoints = Select[blue, Norm[#] < 4.9 &];
initialcond = Table[{q[0] == initialpoints[[j, 1]], p[0] == initialpoints[[j, 2]]}, {j, 1, Length[initialpoints]}];
solutions[equations_] := Table[NDSolve[{equations, initialcond[[j]]}, {q[t], p[t]}, {t, -1, 10}], {j, 1, Length[initialcond]}]
plot[solution_, tmax_, plotlabel_] := Show[
  ParametricPlot[{q[t], p[t]} /. solution, {t, tmax - 0.5, tmax}, PlotStyle -> {Thick},Background -> White,  Axes -> False, PlotRange -> 5.1 {{-1, 1}, {-1, 1}}, PlotLabel -> plotlabel, LabelStyle -> {Black, Bold}, RegionFunction -> Function[{x, y, t}, Sqrt[x^2 + y^2] < 5], ColorFunction -> Function[{x, y, t}, Directive[ColorData["GrayTones"][t/\[Pi]] , Opacity[t^3] ] ]
   ]
  ,
  Graphics[{Black, PointSize[0.02], Point[Select[Flatten[{q[t], p[t]} /. solution /. {t -> tmax}, 1], Norm[#] < 5 &] ], Thick, Circle[{0, 0}, 5]}]
  ]
(*Solve the equations*)
solhyperbolic = solutions[eqs[DiagonalMatrix[{-1, 2}]]];
solelliptic = solutions[eqs[RotationMatrix[\[Pi]/2]]];
solspiralstable = solutions[eqs[-3 RotationMatrix[\[Pi]/5]]];
solspiralunstable = solutions[eqs[1.5*RotationMatrix[\[Pi]/5]]];
(*Plot and animate*)
frames = Table[
   GraphicsGrid[{{
      plot[solspiralstable, \[Tau], "Stable fixed point"], 
      plot[solspiralunstable, \[Tau], "Unstable fixed point"]
      }, {
      plot[solhyperbolic, \[Tau], "Hyperbolic fixed point"], 
      plot[solelliptic, \[Tau], "Elliptic fixed point"]
      }}]
   , {\[Tau], 10^-3, 2, 0.05}];
ListAnimate[frames]