(* Orbits of solutions of {x'[t] == y[t], y'[t] == -x[t] - y[t]/10}
   for a range of initial values. The orbits are color-coded from 
   blue to red as t increases. This example is on p. 67 of "VisualDSolve" 
   by Dan Schwalbe and Stan Wagon  *)

Needs["VisualDSolve`"];

PhasePlot[{x'[t] == y[t], y'[t] == -x[t] - y[t]/10},
  {x[t], y[t]}, {t, 0, 15 Pi}, {x, -1.9, 2.3}, {y, -2, 2},
  InitialValues -> Table[{1, 1 + t}, {t, 0, 1, 0.1}],
  ParametricPlotFunction -> ColorParametricPlot,
  PlotStyle -> AbsoluteThickness[2], AspectRatio -> Automatic,
  WindowShade -> GrayLevel[0.5], DirectionArrow -> False]