Instructor: José Antonio Agapito Ruiz
UCSC Summer Session II, 2003
Math 24

Ordinary Differential Equations

Final Solution

This report does not contain full answers to all questions of the final. Instead, you will find the final is used as an excuse to get more familiar with Mathematica, version 4.2. Not only some relevant Mathematica code is included but also the use of notebooks and some of its nice features is shown. At the same time, I use Problem 1 and 2 to clarify some concepts to you. I hope you will find this useful.

Problem 1.

(a) Find the general solution of the given system of equations and describe the behavior of the solution as t → ∞. Also draw a direction field and plot a few trajectories of the system. Remember to justify your work.

(i) (Overscript[x, ⇀])'=(({{1, -2}, {3, -4}}))(Overscript[x, ⇀]);

Clear[a]
a=(({{1, -2}, {3, -4}}));
len = 
        Length[a];
Solve[
      Det[a - λ IdentityMatrix[len]] == 
        0, λ]

Eigenvalues[a]

Eigensystem[a]

CharacteristicPolynomial[a, λ]

eq = {(Subscript[x, 1])'[t] == 
            (Subscript[x, 1])[t] - 2 (Subscript[x, 2])[t], (Subscript[x, 2])'[t] == 
            3 (Subscript[x, 1])[t] - 4 (Subscript[x, 2])[t]};
sol = DSolve[eq, {(Subscript[x, 1])[t], (Subscript[x, 2])[t]}, t]

It is very instructive to compare this answer with the one we can obtain right away from simple inspection by just looking at the eigenvalues and eigenvectors of the coefficient matrix a. The straightforward answer is

(Overscript[x, ⇀])(t)=(Subscript[C, 1])(e^(-2 t))(({{2}, {3}}))+(Subscript[C, 2])(e^(-t))(({{1}, {1}})),

where (Subscript[C, 1])and (Subscript[C, 2])are arbitrary constants. You may want to apply the commands Simplify or FullSimplify to the answers given by Mathematica for (Subscript[x, 1])and (Subscript[x, 2]) to get simpler expressions. I did that but I got the same answers as before

FullSimplify[sol]

So, (since I don't know at the moment any other way to do it with Mathematica), I performed the simplification by hand and obtained

(Overscript[x, ⇀])(t)=((Subscript[C, 2]) - (Subscript[C, 1]))(e^(-2 t))(({{2}, {3}}))+(3  (Subscript[C, 1]) - 3  (Subscript[C, 2]))(e^(-t))(({{1}, {1}})),

which is equivalent to the expression we obtained in (1). Remember, (Subscript[C, 1]) and (Subscript[C, 2]) are arbitrary constants.

Finally, we draw a direction field and plot a few trajectories of the system for different initial conditions.

Clear[a, b, eq, sol]
eq = {(Subscript[x, 1])'[t] == 
            (Subscript[x, 1])[t] - 2 (Subscript[x, 2])[t], (Subscript[x, 2])'[t] == 
            3 (Subscript[x, 1])[t] - 4 (Subscript[x, 2])[t]};
sol = 
      DSolve[Join[eq, {(Subscript[x, 1])[0] == a, (Subscript[x, 2])[0] == b}], {(Subscript[x, 1])[t], 
            (Subscript[x, 2])[t]}, t] // FullSimplify

We let a vary from -12 to 12 in steps of 2, and then we give b first the value 30 and then the value -30:

Clear[solset]
solset = 
        Table[sol[[1]], {a, -12, 12, 
            2}];
Map[
        ParametricPlot[
            Evaluate[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. solset /. b -> #], {t, 0, 
              6}, PlotStyle -> AbsoluteThickness[0.25], 
            DisplayFunction -> 
              Identity] &, {30, -30}];
Show[%, DisplayFunction -> $DisplayFunction];

These trajectories are not spiral. We can also plot a direction field:

Clear[sol]
sol = 
        DSolve[Join[
            eq, {(Subscript[x, 1])[0] == -10, (Subscript[x, 2])[0] == 20}], {(Subscript[x, 1])[t], 
            (Subscript[x, 2])[t]}, t];
Needs["Graphics`PlotField`"]
Block[{$DisplayFunction = Identity}, 
        g1 = PlotVectorField[{u - 2  v, 3  u - 4  v}, {u, -20, 
              20}, {v, -20, 20}, 
            PlotPoints -> 11]; 
g2 = 
          ParametricPlot[
            Evaluate[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. sol], {t, 0, 
              6}]];
\(\(Show[g1, g2, \ Axes -> True, 
        AxesLabel -> {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])""}];\)\)

(ii)  (Overscript[x, ⇀])'=(({{1, 2}, {-5, -1}}))(Overscript[x, ⇀]);

Clear[a]
a=(({{1, 2}, {-5, -1}}));
len = 
        Length[a];
Solve[
      Det[a - λ IdentityMatrix[len]] == 
        0, λ]

Eigensystem[a]

The general solution is

Clear[eq, sol]
eq = {(Subscript[x, 1])'[t] == 
            (Subscript[x, 1])[t] + 2 (Subscript[x, 2])[t], (Subscript[x, 2])'[t] == -5 (Subscript[x, 1])[t] - 
              (Subscript[x, 2])[t]};
sol = DSolve[eq, {(Subscript[x, 1])[t], (Subscript[x, 2])[t]}, t]

These formulae can be written as

(Overscript[x, ⇀])(t)=(Subscript[C, 1])(e^(-3 i t))(({{-1 + 3  i}, {5}}))+(Subscript[C, 2])(e^(3  i t))(({{-1 - 3  i}, {5}})),

where (Subscript[C, 1]) and (Subscript[C, 2]) are arbitrary constants. We can also write the general solution of the system (ii) in terms of a different set of fundamental solutions, by simply taking the real and imaginary parts of the corresponding fundamental solutions shown in (2). Namely,

FullSimplify[sol]

Recall, if the matrix of coefficients were (({{α, β}, {-β, α}})) then the general solution is of the form

(Overscript[x, ⇀])(t)=(Subscript[C, 1])(e^(α t))(({{Cosβ t}, {-Sinβ t}}))+(Subscript[C, 2])(e^(α t))(({{Sinβ t}, {Cosβ t}})),

where α and β are the real and imaginary parts of the complex eigenvalue λ = α + β i. The nice feature about formula (3) is that we can easily see that its corresponding phase portrait consists of circular trajectories if α = 0 or spirals otherwise. In our case, we don't have this standard matrix , but we can say that the phase trajectories of system (ii) are closed curves already since we do have α = 0. In fact, we can write the expressions for (Subscript[x, 1])and (Subscript[x, 2]) obtained above with Mathematica as

(Overscript[x, ⇀]) (t)=(({{Cos3t + (1/3) Sin3t, (2/3) Sin3t}, {-(5/3) Sin3t, Cos3t - (1/3) Sin3t}}))(({{(Subscript[C, 1])}, {(Subscript[C, 2])}})),

(Overscript[x, ⇀])(t)=(Subscript[C, 1])(({{Cos3t + (1/3) Sin3t}, {-(5/3) Sin3t}}))+(Subscript[C, 2])(({{(2/3) Sin3t}, {Cos3t - (1/3) Sin3t}})).

One can show from here that phase trajectories of system (ii) are ellipses.

Now, we draw a direction field and plot a few trajectories of the system for different initial conditions.

Clear[a, b, eq, sol]
eq = {(Subscript[x, 1])'[t] == 
            (Subscript[x, 1])[t] + 2 (Subscript[x, 2])[t], (Subscript[x, 2])'[t] == -5 (Subscript[x, 1])[t] - 
              (Subscript[x, 2])[t]};
sol = 
      DSolve[Join[eq, {(Subscript[x, 1])[0] == a, (Subscript[x, 2])[0] == b}], {(Subscript[x, 1])[t], 
            (Subscript[x, 2])[t]}, t] // FullSimplify

We let a vary from -12 to 12 in steps of 2, and then we give b first the value 30 and then the value -30:

Clear[solset]
solset = 
        Table[sol[[1]], {a, -12, 12, 
            2}];
Map[
        ParametricPlot[
            Evaluate[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. solset /. b -> #], {t, 0, 
              6}, PlotStyle -> AbsoluteThickness[0.25], 
            DisplayFunction -> 
              Identity] &, {30, -30}];
Show[%, DisplayFunction -> $DisplayFunction];

A plot of a phase trajectory and the direction field can be obtained as follows

Clear[sol]
sol = 
        DSolve[Join[
            eq, {(Subscript[x, 1])[0] == -10, (Subscript[x, 2])[0] == 20}], {(Subscript[x, 1])[t], 
            (Subscript[x, 2])[t]}, t];
Needs["Graphics`PlotField`"]
Block[{$DisplayFunction = Identity}, 
        g1 = PlotVectorField[{u + 2  v, -5 u - v}, {u, -20, 
              20}, {v, -20, 20}, 
            PlotPoints -> 11]; 
g2 = 
          ParametricPlot[
            Evaluate[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. sol], {t, 0, 
              6}]];
\(\(Show[g1, g2, \ Axes -> True, 
        AxesLabel -> {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])""}];\)\)

(iii)  (Overscript[x, ⇀])'=(({{-3, (5/2)}, {-(5/2), 2}}))(Overscript[x, ⇀]);

Clear[a]
a=(({{-3, (5/2)}, {-(5/2), 2}}));
len = 
        Length[a];
Solve[
      Det[a - λ IdentityMatrix[len]] == 
        0, λ]

Eigensystem[a]

We know that (0, 0) is not an eigenvector. This tells us that we need a generalized eigenvector to be able to write down the general solution. First, let's see what Mathematica gives us.

Clear[eq, sol]
eq = {(Subscript[x, 1])'[t] == -3 (Subscript[x, 1])[t] + (5/2) (Subscript[x, 2])[t], (Subscript[x, 2])'[
              t] == -(5/2) (Subscript[x, 1])[t] + 
              2  (Subscript[x, 2])[t]};
sol = DSolve[eq, {(Subscript[x, 1])[t], (Subscript[x, 2])[t]}, t]

FullSimplify[sol]

We can write this as

(Overscript[x, ⇀]) (t)=(Subscript[C, 1])(e^(-(t/2)))(({{1}, {0}}))+((5/2) (Subscript[C, 2]) - (5/2) 
                    (Subscript[C, 1]))(te^(-(t/2)))(({{1}, {1}}))+(Subscript[C, 2])(e^(-(t/2)))(({{0}, {1}})).

If we let (Subscript[C, 2]) = 2  (Subscript[C, 1])and then set (Subscript[K, 1]) ≡ (Subscript[C, 1])and (Subscript[K, 2]) ≡ (5/2) (Subscript[C, 1]), we get

(Overscript[x, ⇀])(t)=(Subscript[K, 1])(e^(-(t/2)))(({{1}, {1}}))+(Subscript[K, 2])[(te^(-(t/2)))(({{1}, {1}}))+(e^(-(t/2)))(({{0}, {(2/5)}}))].

This is the form of the general solution we can get by hand, where (0, (2/5)) is the generalized eigenvector obtained as a solution to the equation

(A - (Subscript[λ, 1]) I) (Subscript[(Overscript[v, ⇀]), 2]) = 
        (Subscript[(Overscript[v, ⇀]), 1]).

In our case, (Subscript[(Overscript[v, ⇀]), 1]) = (1, 1) is the eigenvector corresponding to the eigenvalue (Subscript[λ, 1]) = -(1/2)and (Subscript[(Overscript[v, ⇀]), 2]) is a generalized eigenvector.

Now, we draw phase trajectories and the direction field corresponding to this system.

Clear[a, b, eq, sol]
eq = {(Subscript[x, 1])'[t] == -3 (Subscript[x, 1])[t] + (5/2) (Subscript[x, 2])[t], (Subscript[x, 2])'[
              t] == -(5/2) (Subscript[x, 1])[t] + 
              2  (Subscript[x, 2])[t]};
sol = 
      DSolve[Join[eq, {(Subscript[x, 1])[0] == a, (Subscript[x, 2])[0] == b}], {(Subscript[x, 1])[t], 
          (Subscript[x, 2])[t]}, t]; TableForm[FullSimplify[sol]]

We let a vary from -6 to 6 in steps of 3, and then we give b first the value 10 and then the value -10:

Clear[solset]
solset = 
        Table[sol[[1]], {a, -6, 6, 3}];
Map[
        ParametricPlot[
            Evaluate[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. solset /. b -> #], {t, 0, 
              6}, PlotStyle -> AbsoluteThickness[0.25], 
            DisplayFunction -> 
              Identity] &, {10, -10}];
Show[%, DisplayFunction -> $DisplayFunction];

And a plot for a phase trajectory and the direction field is

Clear[sol]
sol = 
        DSolve[Join[
            eq, {(Subscript[x, 1])[0] == -1, (Subscript[x, 2])[0] == 5}], {(Subscript[x, 1])[t], 
            (Subscript[x, 2])[t]}, t];
Needs["Graphics`PlotField`"]
Block[{$DisplayFunction = Identity}, 
        g1 = PlotVectorField[{-3 u + (5/2) v, -(5/2) u + 
                2  v}, {u, -15, 15}, {v, -15, 15}, 
            PlotPoints -> 11]; 
g2 = 
          ParametricPlot[
            Evaluate[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. sol], {t, 0, 
              6}]];
\(\(Show[g1, g2, \ Axes -> True, 
        AxesLabel -> {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])""}];\)\)

(b) Solve the given value problem and describe the behavior of the solution as t → ∞.

(Overscript[x, ⇀])'=(({{1, 1, 2}, {0, 2, 2}, {-1, 1, 3}}))(Overscript[x, ⇀]),      (Overscript[x, ⇀])(0)=(({{2}, {0}, {1}}))

Clear[a]
a=(({{1, 1, 2}, {0, 2, 2}, {-1, 1, 3}}));
Eigensystem[
      a]

Clear[eq, sol]
eq = {(Subscript[x, 1])'[t] == 
            (Subscript[x, 1])[t] + (Subscript[x, 2])[t] + 2  (Subscript[x, 3])[t], (Subscript[x, 2])'[t] ==  
            2  (Subscript[x, 2])[t] + 2  (Subscript[x, 3])[t], (Subscript[x, 3])'[t] == -(Subscript[x, 1])[t] + 
              (Subscript[x, 2])[t] + 3  (Subscript[x, 3])[t]};
sol = 
      DSolve[Join[eq, {(Subscript[x, 1])[0] == 2, (Subscript[x, 2])[0] == 0, (Subscript[x, 3])[0] == 1}], {(Subscript[x, 1])[t], 
          (Subscript[x, 2])[t], (Subscript[x, 3])[t]}, t]; TableForm[FullSimplify[sol]]

As t → ∞, the solution (Overscript[x, ⇀])(t) goes to ∞ as well.

(c) Find the general solution of the given system of equations.

(Subscript[x, 1])' = -(5/4) (Subscript[x, 1]) + (3/4) (Subscript[x, 2]) + 
        2  t
(Subscript[x, 2])' = (3/4) (Subscript[x, 1]) - (5/4) (Subscript[x, 2]) + (e^t)

Clear[a]
a=(({{-(5/4), (3/4)}, {(3/4), -(5/4)}}));Eigensystem[a]

Clear[eq, sol]
eq = {(Subscript[x, 1])'[
              t] == -(5/4) (Subscript[x, 1])[t] + (3/4) (Subscript[x, 2])[t] + 
              2  t, (Subscript[x, 2])'[
              t] ==  (3/4) (Subscript[x, 1])[t] - (5/4) (Subscript[x, 2])[t] + 
              Exp[t]};
sol = DSolve[eq, {(Subscript[x, 1])[t], (Subscript[x, 2])[t]}, t] // FullSimplify

This can be written as

(Overscript[x, ⇀])(t)=(1/2)((Subscript[C, 1]) + (Subscript[C, 2]))(e^(-(t/2)))(({{1}, {1}}))+(1/2)((Subscript[C, 1]) - (Subscript[C, 2]))(e^(-2 t))(({{1}, {-1}}))+t(({{(5/2)}, {(3/2)}}))-(({{(17/4)}, {(15/4)}}))+(e^t)(({{(1/6)}, {(1/2)}})),

which coincides with the form obtained by hand using the variation of parameters method when setting (Subscript[K, 1]) ≡ (1/2) ((Subscript[C, 1]) + (Subscript[C, 2])) and (Subscript[K, 2]) ≡ (1/2) ((Subscript[C, 1]) - (Subscript[C, 2])).

Problem 2. (Predator-Prey model)

We will study a system  (Overscript[x, ⇀])' = (Overscript[v, ⇀])(
          (Overscript[x, ⇀])), where (Overscript[v, ⇀]) is a nonlinear vector field. The particular system we are interested in is

({{(dx/dt), =, (9 - α x - 3  y) x}, {(dy/dt), =, (-2 + x) y}})

where α ≥ 0 is a parameter. The variable x is the population of prey (rabbits) and y is the population of predators (foxes). We want to understand what happens to these populations as t -> ∞. Our job is to investigate the phase portrait of this system for various values of α in the interval 0 ≤ α ≤ 5. We will follow the steps described in class (see also a nonlinear system analysis, written in Mathematica by Gavin LaRose from the University of Michigan).

Step 1 (Equilibrium points)

We find the equilibrium points by setting both the first and the second components of the vector field (Overscript[v, ⇀]) = ((9 - α x - 3  y) 
            x, (-2 + x) y) equal to zero.

Clear[α, x, y]
Solve[{(9 - α x - 3  y) x == 
          0, (-2 + x) y == 0}, {x, y}]

Clearly, the second equilibrium point above makes sense when α ≠ 0. For α = 0, we get

Clear[x, y]
Solve[{(9 - 3  y) x == 0, (-2 + x) y == 0}, {x, 
        y}]

Step 2 (Linearization. Analysis of eigenvalues for different values of α)

We want to know what type of equilibrium points they are (hyperbolic or non-hyperbolic) If hyperbolic, we know that we get a good approximation of the phase portrait of (6) by studying the linear system (Overscript[y, ⇀])' = 
        D (Overscript[v, ⇀])(p) (Overscript[y, ⇀]) at each hyperbolic equilibrium point p. Thus, we can say that p is like a node, source, center, saddle, and so on, according to the type of equilibrium (0, 0) is in the linear system  (Overscript[y, ⇀])' = 
        D (Overscript[v, ⇀])(p) (Overscript[y, ⇀]). We compute the Jacobian of the vector field (Overscript[v, ⇀]) to then evaluate it at the equilibrium points. Here is a set of commands in Mathematica to perform this task

Clear[α, x, y]
(Subscript[v, 1]) = (9 - α x[t] - 3  y[t]) x[t];  
    (Subscript[v, 2]) = (-2 + x[t]) y[t];
jacobi[f_List, x_List] := 
        Outer[D, f, x];
jac = 
        jacobi[{(Subscript[v, 1]), (Subscript[v, 2])}, {x[t], y[t]}];
MatrixForm[Simplify[jac]]

Again, the equilibrium points are

equi=Solve[{(Subscript[v, 1]) == 0, (Subscript[v, 2]) == 0}, {x[t], 
          y[t]}]

and the eigenvalues of their corresponding Jacobians are

eig=Map[Eigenvalues[jac /. #] &, equi]

Right away we can see that independently of α, the equilibrium point (0, 0) is like a saddle. For the other two equilibrium points, we can try several nonzero values for α to see what we obtain .

α = 1; equi

eig

Looking at their eigenvalues, we see that all three equilibrium points are hyperbolic: (9, 0) is a like a saddle and (2, (7/3)) is like a spiral with converging trajectories.

α = 2; equi

eig

All three equilibrium points are hyperbolic: ((9/2), 0) is like a saddle and (2, (5/3)) is like a spiral with converging trajectories.

α = 3; equi

eig

All three equilibrium points are hyperbolic: (3, 0) is like a saddle and (2, 1) is like a node.

α = 4; equi

eig

All three equilibrium points are hyperbolic: ((9/4), 0) is like a saddle and (2, (1/3)) is like a node.

α = 4.5; equi

eig

Here, we have only two equilibrium points: (0, 0) and (2, 0). And (2, 0) is not hyperbolic

α = 5; equi

eig

Again, all three equilibrium points are hyperbolic: ((9/5), 0) is like a node and (2, -(1/3)) is like a saddle.

When α = 0, we can independently compute eigenvalues for the corresponding equilibrium points. Now, I denote (Subscript[v, 1, 0]) and (Subscript[v, 2, 0]) the first and second components of the vector field (Overscript[v, ⇀]).

Clear[x, y, equi]
(Subscript[v, 1, 0]) = (9 - 3  y[t]) x[t];  
    (Subscript[v, 2, 0]) = (-2 + x[t]) y[t];
jac = 
        jacobi[{(Subscript[v, 1, 0]), (Subscript[v, 2, 0])}, {x[t], 
            y[t]}];
MatrixForm[Simplify[jac]]

equi=Solve[{(Subscript[v, 1, 0]) == 0,  
          (Subscript[v, 2, 0]) == 0}, {x[t], y[t]}]

eig=Map[Eigenvalues[jac /. #] &, equi]

Thus, for α = 0, we have two equilibrium points: (0, 0) is hyperbolic, whereas (2, 3) is not since its eigenvalues have real part equal to zero.

Step 3 (Direction fields, phase portraits)

First, we plot direction fields of the system (6) for values of α = 0, 1, 2, 3, 4, 5, to get started.

Clear[α]
<< Graphics`PlotField`
Table[
        PlotVectorField[{(9 - α x - 3  y) x, (-2 + x) 
              y}, {x, -1, 5}, {y, -1, 5}, Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          PlotLabel -> SequenceForm["α=", α], 
          PlotPoints -> 20], {α, 0, 
          5}];

These plots give us already a good picture of what is going on with sytem (6) for different values of the parameter α. By simple inspection we can guess where the equilibrium points are and the type of them (something we have done already).

It is always a good idea to plot the curves where x' and y' vanish (some authors call these curves nullclines). For instance, if you are doing things by hand, you can proceed by steps to construct a good sketch of the direction field of (6). For x' = 0, the resulting curves determine regions on the phase plane. In each region you analyze the horizontal component of the vector field (Overscript[v,  ⇀])  in (6). Along these curves, the vector field (Overscript[v,  ⇀])  is vertical since the horizontal component of (Overscript[v,  ⇀])  is zero (x' = 0). For y' = 0, you do a similar analysis but now focusing on the vertical components of (Overscript[v,  ⇀]) . Obviously, here, the vector field (Overscript[v,  ⇀])  is horizontal along the corresponding curves where y' = 0. Finally, you overlap the two plots above to get a sketch of the direction field of the system. We do this for α = 0, 1, 2, 3, 4, 5. Red curves represent points where x' = 0 and blue curves are points where y' = 0.

Clear[α]
h1=y[t] /. Solve[(Subscript[v, 1]) == 0, y[t]][[1]]/.x[t]->x

v2 = y[t] /. Solve[(Subscript[v, 2]) == 0, y[t]][[1]] /. x[
            t]->x

\(\(Table[
        Block[{$DisplayFunction = Identity}, 
          hor1 = Plot[h1, {x, 0, 10}, PlotStyle -> RGBColor[1, 0, 0]]; 
          hor2 = Graphics[{RGBColor[1, 0, 0], Line[{{0, \(-5\)}, {0, 5}}]}]; 
          ver1 = Graphics[{RGBColor[0, 0, 1], Line[{{2, \(-5\)}, {2, 5}}]}]; 
          ver2 = Plot[v2, {x, 0, 10}, PlotStyle -> RGBColor[0, 0, 1]]; 
          equipts = 
            Graphics[{PointSize[ .02], Point[{0, 0}], 
                Point[{9\/α, 0}], 
                Point[{2, \(9 - 2  α\)\/3}]}]]; 
Show[hor1, hor2, ver1, ver2, equipts, Axes -> True, 
          PlotLabel -> SequenceForm["α=", α], 
          AxesLabel -> {"x", "y"}], {α, 1, 
          5}];\)\)

The intersections of red and blue curves are the equilibrium points we have previously obtained. Now, we plot some phase trajectories for different initial conditions.

\(Clear[α, eq, sol]
 (*\ list\ of\ equations\ *) \)
\(eq = {\(x'\)[t] == v\_1, \(y'\)[t] == v\_2}; \ 
 (*\ 
      10\ initial\ conditions\ *) \ \)
\(iclist = {{0.1, 0.05}, {1, 2}, {2.5, 1}, {3, 0.05}, {3.8, 3}, {4, 
          0}, {4.3, 1}, {4.54, 0.2}, {5, 0.1}, {6, 1}}; 
    Table[Block[{$DisplayFunction = Identity}, 
        getPlot[eq_, x0_, y0_] := 
          Module[{sol, xs, ys}, 
            sol = NDSolve[{eq, x[0] == x0, y[0] == y0}, \ {x[t], 
                  y[t]}, {t, 0, 50}]; {xs, ys} = {x[t], y[t]} /. 
                sol[\([\)\(1\)\(]\)]; 
            ParametricPlot[{xs, ys}, {t, 0, 50}, 
              DisplayFunction -> Identity]]; 

        trajectories = 
          Table[getPlot[eq, iclist[\([\)\(i, 1\)\(]\)], 
              iclist[\([i, 2]\)]], {i, 1, Length[iclist]}]]; 

      trajplot = 
        Show[trajectories, DisplayFunction -> $DisplayFunction, 
          AspectRatio -> 1, PlotRange -> All, 
          PlotLabel -> SequenceForm["α=", α], 
          AxesLabel -> {"x", "y"}], {α, 1, 
        5}]\)

Finally, we can put all the plots together: direction fields, nullclines, equilibrium points and some phase trajectories. I leave this for you as an exercise. You may want to use again the command Table to plot all these graphs for different values of α, all at once.

Step 4 (Other type of diagrams: x vs. t, y vs. t)

Here, just follow the pattern as in the Mathematica tutorial (Part b, pp. 674-679).

Problem 3. (Extra credit)

Solve Problem 18 of Section 7.8 of the textbook.

The solution is similar to that of  Problem 17 in the same section. Just follow the steps.