<< Graphics`PlotField`

example 5.7  p.68

eq = y'[x] == -2 x y[x]

yx1 = 
      y[x] /. DSolve[{eq, y[0] == 1}, y[x], 
            x][[1]]
yx2 = 
      y[x] /. DSolve[{eq, y[0] == 2}, y[x], 
            x][[1]]
yx3 = 
      y[x] /. DSolve[{eq, y[0] == 4}, y[x], 
            x][[1]]
yx4 = 
      y[x] /. DSolve[{eq, y[0] == 8}, y[x], x][[1]]

zu0 = 
      Plot[{yx1, yx2, yx3, yx4}, {x, -3, 3}, 
        PlotRange ->  {0, 10}]

zu1 = 
      PlotVectorField[{1, -2 x y}, {x, -3, 3}, {y, 0, 
          10}]

Show[zu0, zu1]

example 5.8  p.68

eq = h'[t] == -k Sqrt[h[t]] + q

para = {k ->  0.1, h0 ->  1, q ->  0.05}

h1 = 
      h[t] /. NDSolve[{eq, h[0] == h0} /. para, 
          h[t], {t, 0, 100}]
h2 = 
      h[t] /. NDSolve[{eq, h[0] == 0.5} /. para, 
          h[t], {t, 0, 100}]
h3 = 
      h[t] /. NDSolve[{eq, h[0] == 0.1} /. para, 
          h[t], {t, 0, 100}]
h4 = 
      h[t] /. NDSolve[{eq, h[0] == 1.5} /. para, 
          h[t], {t, 0, 100}]

zu0 = 
      Plot[{h1, h2, h3, h4}, {t, 0, 80}, 
        PlotRange ->  {0, 1}]

vec = {1, -k Sqrt[h] + q} /. para

zu2 = PlotVectorField[vec, {t, 1, 80}, {h, 0, 1}]

Show[zu0, zu2]

example 5.9  p.70

eq = v'[r] v[r] == -g R^2/r^2

ans1 = 
        DSolve[{eq, v[R] == v0}, v[r], r];
vr = v[r] /. ans1[[2]]

para = {R ->  6372 10^3, g ->  9.8}

v1 = vr/10^3 /. r ->  R k /. para /. 
        v0 ->  10 10^3

v2 = vr/10^3 /. r ->  R k /. para /. 
        v0 ->  11.2 10^3

zu0 = 
      Plot[{v1, v2}, {k, 1, 10}, PlotRange ->  {0, 15}]

vec = {v/R, -g R^2/r^2/10^3} /. {v ->  y 10^3, 
            r ->  R k} /. para

zu3 = PlotVectorField[vec, {k, 1, 10}, {y, 0, 14}]

Show[zu0, zu3]