(* Phase plot of {x'[t] == x[t] - y[t]^2Cos[y[t]], y'[t] == - y[t] + x[t] Sin[x[t]]}.
   Shaded regions show where the vector field is pointing in a common direction;
   borders between regions are nullclines (lines of zero slope). "Fish" shapes
   indicate the direction and magnitude of the vector field over short
   trajectories, color-coded for magnitude. Black dots show equilibrium points
   (x' = y' = 0). This example is on p. 84 of "VisualDSolve" by Dan Schwalbe and 
   Stan Wagon  *)

Needs["VisualDSolve`"];

PhasePlot[
    {x'[t] == x[t] - y[t]^2Cos[y[t]], y'[t] == - y[t] + x[t] Sin[x[t]]},
    {x[t], y[t]}, {x, -10, 10}, {y, -10, 10},
    Nullclines -> True, NullclineShading -> True,
    ShowEquilibria -> True, 
    EquilibriumPointStyle -> {GrayLevel[0], PointSize[0.01]},
    NullclinePlotPoints -> 80, PlotPoints -> 500,
    FlowField -> True, FieldLogScale -> Automatic,
    FieldLength -> 1.0, FieldThickness -> 0.6, FieldColor -> Hue];