DE.05 First Order Differential Equations
 BASICS

Needs["Graphics`Colors`"];
Needs["Graphics`FilledPlot`"];

If[MemberQ[{"Global`","Graphics`Common`GraphicsCommon`","Vector3D`"}, 
    Context[ScaleFactor]],Remove["*`ScaleFactor"]];
<<"Graphics`Legend`";
<<"Graphics`Common`GraphicsCommon`";
Remove["*`ScaleFactor"];
Graphics`Colors`GosiaGreen=RGBColor[0, 0.392187, 0];
Off[General::spell];
Off[General::spell1];
Off[ParametricPlot3D::ppcom];
Off[ParametricPlot::ppcom];
Off[Plot::plnr];
Needs["Calculus`Master`"];
Off[Unique::"usym"];
SetOptions[Limit,Analytic->True];

CMView = {2.7, 1.6, 1.2};

Off[Graphics::gpt]

New Vector Drawers

(* :Title:  NewVectorDraw  *)
(* :Author: Bruce Carpenter *)
(* :Summary:  Construct 2 and 3 dimensional vectors,
     arrows, and arrowheads as graphics 
     objects.
*)
(* :Date:  Copyright 1997, Math Everywhere, Inc. 
*)
(*  :Drawbacks:  Since the vectors are drawn out of 
     context, the shapes of the arrowheads are only 
     correct with AspectRatio set to Automatic.
*)
BeginPackage["Vector3D`","Graphics`Colors`"];

Axes3D::usage = "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8].";

Perpend::usage = "Perpend[a] returns a unit vector which is perpendicular to the vector a.";

Arrow::usage = "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. Arrow[a, Tail -> tail] gives a vector from tail to a.";
	
Vector::usage = "Vector[a, b] gives a 2D or 3D vector from a to b.";

ArrowHead::usage = "ArrowHead[a, vec] produces an arrowhead with its tip placed at the point a, pointing in the direction of the vector vec.";

Tail::usage = "Tail -> point puts the tail of the vector at point.";

Norm::usage = "Norm[vec] returns the Euclidean length of vec.";

Normalize::usage = "Normalize[vec] returns a unit vector in the direction of vec, provided vec is nonzero.";

Aperture::usage = "Aperture is the ratio of the radius of the base to the length of the head of a vector.";

TipSize::usage = "TipSize specifies an absolute size for the inner tip length of the head of a vector.";

TipRatio::usage = "TipRatio is the ratio of the inner length of the arrowhead to the length of the shaft of a vector.";

HeadRatio::usage = "HeadRatio is the ratio of the outer length of the arrowhead to the length of the shaft of a vector.";

HeadSize::usage = "HeadSize specifies an absolute size for the head of a vector.";

ScaleFactor::usage = "ScaleFactor specifies the amount to scale a vector in length and may be either a number or a pure function.";

ZeroVectorPointSize::usage = "ZeroVectorPointSize is the size of the point used to represent a zero vector.";

Options[Vector] = Options[Arrow] = Options[ArrowHead] =
{HeadRatio -> 0.18, TipRatio -> 0.14, Aperture -> 0.3, ScaleFactor -> 1, HeadType -> Polygon, ShaftQ -> True, EdgesQ -> True, VectorColor -> RGBColor[0,0,1], ShaftWidth -> 0.005, ZeroVectorPointSize -> 0.01};
		
Begin["`Private`"];

SetOptions[ParametricPlot,AspectRatio -> Automatic];
SetOptions[Plot,AspectRatio -> Automatic];
SetOptions[Graphics,AspectRatio -> Automatic];

Axes3D[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u/3,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u/3,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,-u/3},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
Axes3D[u_] := Axes3D[u,u/8];
				
Norm[a_?VectorQ] := N[Sqrt[a.a]];

Normalize[a_?VectorQ] := If[a.a == 0, a, a/Norm[a]];

Perpend[a_?VectorQ] := Normalize[{a[[2]],-a[[1]]}]/;Length[a]==2
Perpend[a_?VectorQ] := Normalize[If[a[[1]]==0,{1,0,0},{0,a[[3]],-a[[2]]}]]/;Length[a]==3

Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/8.]},{k,0,8}]; 

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, tip},
   {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, TipRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]; tr = 0.75 hr];
   If[NumberQ[ts], tr = ts/Norm[diff]];
   Graphics[{If[sq,{vc,Thickness[sw],Line[{from,tip - tr diff}]},{}],
      With[{trans = ap hr Norm[diff] Perpend[diff]},
      {vc, Thickness[sw], ht[{tip,
                   (tip - hr diff + trans),
                   (tip - tr diff),
                   (tip - hr diff - trans),tip}]}]}]
] /; Length[from]==Length[to]==2

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, eq, sw, hs, hr, tr,ap, sf, vc, zvps, diff = to - from, tip},
   {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics3D[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]];
   tr = hr/2;
   Graphics3D[{If[sq,{vc,Thickness[sw],Line[{from,tip-tr diff}]},{}], If[eq, {}, EdgeForm[]], {SurfaceColor[vc],
      With[{temp = Perpend[diff]},
      Polygon[Append[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[a_, Tail->b_, opts___] := Vector[b,a+b,opts]

Arrow[a_, opts___] := Vector[Table[0,{Length[a]}],a,opts]

ArrowHead[a_, b_, opts___] := Vector[a-b,a,ShaftQ->False,opts];

End[];

EndPackage[];
Null

VectorField and Phase Lines

(*  :Title: VectorAndPhase.m  *)

(*  :Author:  Bill Davis  *)

(*  :Dates:  about March, 1996 *)

(*  :Drawbacks:  I wish someone serious would come and turn these into
                 Really carefully done packages.  *)

BeginPackage["PhaseStuff`","Vector3D`",
		"Utilities`FilterOptions`"];

PhaseLine::usage = "PhaseLine[fcn,{var,a,b},tint,pos]
	makes a vertical phase line for the function, fcn, of the variable, var,
	for the variable going from a to b, with 
	a color given by tint,  at the x-coordinate pos.";
	
PhaseLine3D::usage = "PhaseLIne3D[fcn,{var,a,b},tint,loc]
	makes a vertical phase line for the function, fcn, of the variable, var,
	for the variable going from a to b, with 
	a color given by tint, and base point at loc = {x,y}";
	
VectorField::usage = "VectorField[{function,{var1,a,b,space1},{var2,c,d,space2}]
	will make a vector field of normalized vectors, colored by
	magnitude, for the function of the variables on the
	range {r1,r2} with spacing between points given by 
	{space1, space2}.  If a single function of a single variable is given,
	the vector field {1,f} will be plotted.";
	
StikPikXZ::usage = "StikPikXZ[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
StikPikYZ::usage = "StikPikYZ[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
StikPikXY::usage = "StikPikXY[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
Begin["`Private`"];	

PhaseLine[fcn_,{var_,a_,b_},tint_,pos_:0] := 
	Module[{len,len1,k,zset,zero,zero1,signs,stick,dots,pointers},
	zset = Solve[fcn==0,var];
	zero = If[Length[zset]==0,{},var/.zset];
	zero = Sort[zero];
	zero = Select[zero,((# >= a)&&(# <= b))&];
	zero1 = Sort[Union[zero,{a,b}]];
	len = Length[zero];
	len1 = Length[zero1];
	signs = Table[Sign[fcn/.var->((zero1[[k+1]] + zero1[[k]])/2)],
		{k,1,len1-1}];
	stick = {Thickness[0.007],Line[{{pos,a},{pos,b}}]};
	dots = Point[{pos,#}]&/@zero;
	dots = Join[{PointSize[0.05],tint},dots];
	heads = Table[MakeHead[zero1[[k]],zero1[[k+1]],pos,
		signs[[k]],b-a,tint],{k,1,len1-1}];
	Return[{Graphics[Join[{tint},stick,dots]],heads}]
		]
		
MakeHead[a_,b_,c_,s_,slen_,tint_] := Module[{u,v,w},
	u = Min[2(b-a),3/4 slen];
	Return[ArrowHead[{c,a(.5-s/5)+b(.5+s/5)},
		{0,s u},VectorColor->tint]]]
		
PhaseLine3D[fcn_,{var_,a_,b_},tint_,pos_:{0,0}] := 
	Module[{len,len1,k,zset,zero,zero1,signs,stick,dots,pointers},
	zset = Solve[fcn==0,var];
	zero = If[Length[zset]==0,{},var/.zset];
	zero = Sort[zero];
	zero = Select[zero,((# >= a)&&(# <= b))&];
	zero1 = Sort[Union[zero,{a,b}]];
	len = Length[zero];
	len1 = Length[zero1];
	signs = Table[Sign[fcn/.var->((zero1[[k+1]] + zero1[[k]])/2)],
		{k,1,len1-1}];
	stick = {Thickness[0.007],Line[{{pos[[1]],pos[[2]],a},
		{pos[[1]],pos[[2]],b}}]};
	dots = Point[{pos[[1]],pos[[2]],#}]&/@zero;
	dots = Join[{PointSize[0.03],tint},dots];
	heads = Table[MakeHead3D[zero1[[k]],zero1[[k+1]],pos,
		signs[[k]],b-a,tint],{k,1,len1-1}];
	Return[{Graphics3D[Join[{tint},stick,dots]],heads}]
		]
		
MakeHead3D[a_,b_,c_,s_,slen_,tint_] := Module[{u,v,w},
	u = Min[2(b-a),1/2 slen];
	Return[ArrowHead[{c[[1]],c[[2]],a(.5-s/5)+b(.5+s/5)},
		{0,0,s u},VectorColor->tint]]]
		
		
VectorField[function_List,{var1_,a_,b_,sp1_:1},{var2_,c_,d_,sp2_:1}] := 
	Module[{i,j,U,V,UU,VV,ulow,uhigh,vlow,vhigh,maxval,minval,
		lengths,directions,locations,mlen},
	mlen = Min[sp1,sp2];
	normf[uu_,vv_] = Sqrt[function[[1]]^2+function[[2]]^2]/.
		{var1->uu,var2->vv};
	lengths = Table[N[normf[U,V]],{U,a,b,sp1},{V,c,d,sp2}];
	maxval = Max[lengths]; minval = Min[lengths];
	If[minval == maxval,maxval += 1];
	locations = Table[{uu,vv},{uu,a,b,sp1},
		{vv,c,d,sp2}];
	directions = Table[If[N[normf[U,V]]==0,{0,0},
		mlen (function/.{var1->U,var2->V})/(1.3 N[normf[U,V]])],
			{U,a,b,sp1},{V,c,d,sp2}];
	{UU,VV} = Dimensions[directions][[{1,2}]];
	Return[Table[Arrow[directions[[i,j]],
			Tail->locations[[i,j]],Aperture->1,
			HeadRatio->.4,ShaftWidth->0.008, VectorColor->
			RGBColor[(lengths[[i,j]] - minval)/(maxval-minval),0,
				(maxval -lengths[[i,j]])/(maxval-minval)]],{i,UU},{j,VV}]]
		]
		
VectorField[function_,{var1_,a_,b_,sp1_:1},{var2_,c_,d_,sp2_:1}] :=
	VectorField[{1,function},{var1,a,b,sp1},{var2,c,d,sp2}]
	
	
Needs["Utilities`FilterOptions`"]

FixOptions[a_] := Module[{picture = a, opt3},
If[Options[picture] != {}, 
	opt3 = FilterOptions[Graphics3D,#]&/@Options[picture];
	(* Fix things that might carry 2D characteristics. *)
	opt3 = opt3//. {(AxesLabel -> _) ->(AxesLabel -> None),
			(PlotRange->_) -> (PlotRange->Automatic),
			(PlotRegion -> _) -> (PlotRegion -> Automatic),
			(Ticks -> _) -> (Ticks -> Automatic)};
	picture = Apply[List,picture];
	picture[[-1]] = opt3;
	picture = Apply[Graphics3D,picture]];
	Return[picture]
	]

StikPikXZ[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{u,a,v}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]


StikPikYZ[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{a,u,v}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]


StikPikXY[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{u,v,a}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]
		
End[];

EndPackage[];
Null

ThreeAxes[u,v]

ThreeAxes[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,0},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
ThreeAxes[u_] := ThreeAxes[u,u/8];
ThreeAxes::usage = "ThreeAxes[a,b] makes a standard cartesian axis graphics object
	with x, y, and z running from -a to a, and with axis
	labels b units beyond the tips of the axes.  ThreeAxes[a]
	is ThreeAxes[a,a/8].";

Crisscross

Clear[x];
crisscross = Plot[{Cos[x],Sin[x]},{x,Pi/8,3 Pi/8},
    PlotStyle->{{RGBColor[1,0,0],Thickness[0.02]}},
    Axes->None,AspectRatio->1/3,
    DisplayFunction->Identity];

B.1)  Reading an autonomous diffeq through phase lines

B.1.a.i)

Here's an old friend, the logistic differential equation:

\(Clear[diffeq, y, t, f, starter]; 
r\  = \ 0.45; 
b\  = \ 6; 

    f[t_, y_]\  = \ r\ y\ \((1\  - \ y/b)\); 


    \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

Here is a flow plot shown with plots of four solutions of this diffeq.

{ylow, yhigh} = {0, 9}; {tlow, thigh} = {0, 12}; flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, HeadSize -> 0.6], {t, tlow, thigh, 
          ((thigh - tlow)/8)}, {y, ylow, yhigh, ((yhigh - ylow)/9)}]; 
    \(Clear[y1, y2, y3, y, t]; 
    
{starter1, starter2, starter3, starter4} = {0, 3.0, 6.0, 8.5}; 

    \(endtime = thigh; \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; y3[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y4[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter4, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t], y3[t], y4[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 

    \(setup = 
      Show[solutionplots, flowplot, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)\)

Look at this embellishment of the graphic immediately above.

phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Blue, -1]; 
    flowphaseplot = Show[setup, phaseline];

Explain what those arrowheads on the far left mean.

Take another look:

Show[flowphaseplot];

Shade in to the right of the down arrowhead:

\(downshadow = 
      Graphics[{GrayLevel[0.7], 
          Polygon[{{tlow, 6}, {tlow, yhigh}, {thigh + 1, yhigh}, {thigh + 1, 
                6}}]}]; 


    Show[downshadow, flowphaseplot, Axes -> True, 
      AxesLabel -> {"t", "y"}, PlotRange -> {ylow, yhigh}, 
      AspectRatio -> 1\/GoldenRatio]; \)

The down arrowhead indicates that any solution passing through the shaded part above is going down and must level out at the cutoff point at  y = 6.

Now shade in to the right of the up arrowhead:

\(upshadow = 
      Graphics[{GrayLevel[0.7], 
          Polygon[{{tlow, ylow}, {tlow, 6}, {thigh + 1, 6}, {thigh + 1, 
                ylow}}]}]; 


    Show[upshadow, flowphaseplot, Axes -> True, 
      AxesLabel -> {"t", "y"}, PlotRange -> {ylow, yhigh}, 
      AspectRatio -> 1\/GoldenRatio]; \)

The up arrowhead indicates that any solution passing through the shaded part above is going up and must level out at the cutoff point at  y = 6.

The solution corresponding to y[0] = 6 stays flat. 
All other solutions plot out entirely above or entirely below the flat solution.

B.1.a.ii)

What's the idea behind the phase line?

See how the diffeq was entered:

\(Clear[diffeq, y, t, f, starter]; 
b\  = \ 6; 

    f[t_, y_]\  = \ 0.45\ y\ \((1 - y\/b)\); 


    \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

This allows you to read off y' as a function of y:

yprime = f[t, y]

Plot y'as a function of y:

phaseplot = 
      Plot[f[t, y], {y, ylow, yhigh}, 
        PlotStyle -> {{Thickness[0.01], Blue}}, 
        AxesLabel -> {"y", "y'"}, 
        AspectRatio -> (1/GoldenRatio)];

This plot tells you that 
        y'=f[t, y] > 0 
for y between 0 and 6.0.

The upshot: 
If a solution is ever between 0 and 6.0, then that solution is going up and will level off as it gets near 6.0.

Take another look:

Show[phaseplot];

This plot also tells you that 
        y'=f[t, y] < 0 
for y bigger than 6.0.
The upshot: 
If a solution is ever bigger than 6.0, then that solution is going down and will level off as it gets near 6.0.

This is exactly the same thing that the phase line plot is telling you:

\(\(\(Show[phaseline, Axes -> True, 
      AxesLabel -> {"t", "y"}, 
      PlotRange -> {{tlow - 1.5, thigh}, {ylow, yhigh}}, 
      AspectRatio -> 1\/GoldenRatio]; \)
\)\)

Show[flowphaseplot];

B.1.b.i)

Here's a new diffeq:

\(Clear[diffeq, y, t, f, starter]; 
b\  = \ 6; 

    f[t_, y_]\  = \ 
      0.08\ \((y\  - \ 2.0)\)\ \((4.0\  - \ y)\)\ \((y\  - \ 8.0)\); 


    \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

Here is the phase line:

{ylow, yhigh} = {-1, 10}; {tlow, thigh} = {0, 12}; \(phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Blue, \(-1\)]; 


    phaselineplot = 
      Show[phaseline, PlotRange -> {{tlow - 1.3, thigh}, {ylow, yhigh}}, 
        Axes -> True, AxesLabel -> {"t", "y"}]; \)

Milk that plot for all the information you can get about how solutions behave.

This gives you just about as much information as the flow field.
Take another look:

Show[phaselineplot];

Read off:
Solutions starting with  y[0]  >  8  go down and level off at y = 8.
Solutions starting with  4 < y[0] < 8  go up and level off at y = 8.
Solutions starting with  2 < y[0] < 4  go down and level off at y  = 2
Solutions starting with  y[0] < 2  go up and level off at y = 2.

Check it out with some sample solution plots:

\(Clear[y1, y2, y3, y4, y, t]; 
    
{starter1, starter2, starter3, starter4} = {\(-1.0\), 3.8, 4.2, 10}; 

    \(endtime = thigh; \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; y3[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}][[1
          ]]; y4[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter4, 
            y[t], {t, 0, endtime}][[1
          ]]; \(solutionplots = 
      Plot[{y1[t], y2[t], y3[t], y4[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 


    Show[solutionplots, phaselineplot, Axes -> True, 
      AxesLabel -> {"t", "y"}, PlotRange -> {ylow, yhigh}, 
      AspectRatio -> 1\/GoldenRatio, 
      DisplayFunction -> $DisplayFunction]; \)

Sure enough.
Solutions starting with  y[0] > 8  go down and level off near y = 8.
Solutions starting with  4 < y[0] < 8  go up and level off near y = 8. 
Solutions starting with  2 < y[0] < 4  go down and level off at y = 2.
Solutions starting with   y[0] < 2  go up and level off at y = 2.

B.1.b.ii)

Stay with the same diffeq and look at the phase line plot again:

Show[phaselineplot];

How does the phase plot signal you to beware of extreme sensitivity to errors in starter data for starting values y[0] near y[0] = 4.0?

Take another look:

Show[phaselineplot];

Just above y = 4, you see an up arrowhead and just below y = 4, you see a down arrowhead. 
This tells you that that for y[0] near 4.0, you can expect sensitive dependence on starter data:
In fact, solutions starting at 4.01 go up, but solutions starting at 3.99 go down.
See it happen for two solutions y1[t] and y2[t] with y1[0] = 3.99 and y2[0] = 4.01:

\(Clear[y1, y2, y3, y4, y, t]; 
{starter1, starter2} = {3.99, 4.01}; 

    \(endtime = thigh; \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; \(solutionplots = 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 


    Show[solutionplots, phaselineplot, Axes -> True, 
      AxesLabel -> {"t", "y"}, PlotRange -> {ylow, yhigh}, 
      AspectRatio -> 1\/GoldenRatio, 
      DisplayFunction -> $DisplayFunction]; \)

That's extreme sensitivity to errors in the starting value of y[0] for y[0] near 4.0.

B.2)  For non-autonomous diffeqs, a single phase line is not all that useful; 
         multiple phase lines are sometimes useful

B.2.a.i) Sometimes a single phaseline won't cut it

Here's a diffeq related to the logistic diffeq.

\(Clear[diffeq1, y, t, f1, starter]; 
b\  = \ 7.0; 

    f1[t_, y_]\  = \ 0.2\ y\^2\ \((1 - y\/b)\); 


    \((diffeq1\  = \ {\(y'\)[t]\  == \ f1[t, y[t]], y[0]\  == \ starter})
        \) // ColumnForm\)

Here's a flow plot shown with plots of four solutions of this diffeq.

{ylow, yhigh} = {0, 10}; {tlow, thigh} = {0, 12}; \(flowplot1 = 
      Table[Arrow[{1, f1[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> Normalize, 
          HeadSize -> 0.6], {t, tlow, thigh, \(thigh - tlow\)\/10}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/10}]; 

Clear[y1, y2, y3, y, t]; 
    
{starter1, starter2, starter3, starter4} = {0, 2.0, 5.0, 9.5}; 

    \(endtime = thigh + 1; \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq1 /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; y2[t_] = 
      y[t] /. NDSolve[diffeq1 /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; y3[t_] = 
      y[t] /. NDSolve[diffeq1 /. starter -> starter3, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y4[t_] = 
      y[t] /. \(NDSolve[diffeq1 /. starter -> starter4, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots1 = 
      Plot[{y1[t], y2[t], y3[t], y4[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 

    phaseline1 = PhaseLine[f1[t, y], {y, ylow, yhigh}, Blue, \(-1\)]; 


    setup1 = Show[solutionplots1, flowplot1, phaseline1, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)

The single phase line works beautifully and reveals just about all you need to know about this diffeq. 
Now look at this new diffeq:

\(Clear[newdiffeq, y, t, f2, starter]; 
f2[t_, y_]\  = \ \ y\  - \ \ t; 
    

\((diffeq2\  = \ {\(y'\)[t]\  == \ f2[t, y[t]], y[0]\  == \ starter})
        \) // ColumnForm\)

And look at this flow plot shown with plots of three solutions of this diffeq:

{ylow, yhigh} = {0, 9}; {tlow, thigh} = {0, 6}; flowplot2 = 
      Table[Arrow[{1, f2[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> Normalize, 
          HeadSize -> 0.4], {t, tlow, thigh, ((thigh - tlow)/10)}, {y, 
          ylow, yhigh, ((yhigh - ylow)/10)}]; \(Clear[y1, y2, y3, y, t]; 
    
{starter1, starter2, starter3} = {0.4, 0.99, 3.0}; 
\(endtime = thigh; 
    \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq2 /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; y2[t_] = 
      y[t] /. NDSolve[diffeq2 /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y3[t_] = 
      y[t] /. \(NDSolve[diffeq2 /. starter -> starter3, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots2 = 
      Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 


    setup2 = Show[solutionplots2, flowplot2, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)

Explain why a single phase line wouldn't have been useful for the second diffeq.

Look at both plots side by side:

both = Show[GraphicsArray[{setup1, setup2}]];

On the left, the single phase line does a great job. 
But on the right a single phase line won't work.
Reason: On the right, some of the solutions go up for a while and then turn around and go down. This makes an overall single phase line out of the question.

B.2.a.ii)

Why did that happen?

Look at the way the two diffeqs were entered:

\(Clear[diffeq1, y, t, f1, starter]; 
b\  = \ 7.0; 

    f1[t_, y_]\  = \ 0.2\ y\^2\ \((1 - y\/b)\); 


    \((diffeq1\  = \ {\(y'\)[t]\  == \ f1[t, y[t]], y[0]\  == \ starter})
        \) // ColumnForm\)

\(Clear[diffeq2, y, t, f2, starter]; 
f2[t_, y_]\  = \ \ y\  - \ \ t; 

    
\((diffeq2\  = \ {\(y'\)[t]\  == \ f2[t, y[t]], y[0]\  == \ starter})
        \) // ColumnForm\)

Now look at the vectors that define each flow:

flow1 = {1, f1[t, y]}

flow2 = {1, f2[t, y]}

The slope of the flow vectors for diffeq1 depends only on y.
The slope of the flow vectors for diffeq2 depends on BOTH t and y. See it:

Show[both];

That's why a single phase line is not all that useful for diffeq2.

B.2.a.iii) Recognizing autonomous diffeqs

Folks say that a diffeq is autonomous if a single phase line tells all.
How do you recognize in advance whether a given diffeq is autonomous?

Look at the function f[t, y] used to enter the diffeq. 
If f[t, y] has no dependence on t, the you have an autonomous diffeq. 
Otherwise you don't.
Samples:

\(Clear[diffeq, b, a, y, t, f, starter]; 


    f[t_, y_]\  = \ a\ y\^3\ \((1 - y\^2\/b)\); 


    \((diffeq1\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

Look at:

f[t, y]

No dependence on t. 
The upshot: This diffeq is autonomous and is ripe for analysis through a phase line.

Another:

\(Clear[diffeq, b, a, y, t, f, starter]; 


    f[t_, y_]\  = a\ y\^2\ \((Sin[t] - y\/b)\); 

    \((diffeq1\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

Look at:

f[t, y]

Vivid dependence on t. 
The upshot: This diffeq is not autonomous.

B.2.a.iv) Multiple phaselines for non-autonomous diffeqs

Go with this new non-autonomous diffeq:

\(Clear[newdiffeq, y, t, f, starter]; 
f[t_, y_]\  = y - 0.5\ t\^2; 


    \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

Look at this flow plot shown with plots of three solutions of this diffeq:

{ylow, yhigh} = {0, 12}; {tlow, thigh} = {0, 6}; scalefactor = 0.6; flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.5], {t, tlow, thigh, ((thigh - tlow)/8)}, {y, 
          ylow, yhigh, ((yhigh - ylow)/9)}]; \(Clear[y1, y2, y3, y, t]; 
    
{starter1, starter2, starter3} = {0.4, 0.9, 2}; 
\(endtime = thigh; 

    
y1[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 

    \(y2[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; \)\)\)\(y3[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 


    setup = Show[solutionplots, flowplot, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)

This diffeq is not autonomous. Thus, it is not possible to use a single phase line to analyze it. 
Can you use multiple phase lines to any profit?

Ya sure! You bettcha! 
Make a different phase line for each t and then plot it at the corresponding t like this:

\(Clear[phaseline]; 

    phaseline[t_] := PhaseLine[f[t, y], {y, ylow, yhigh}, GrayLevel[0.5], t]; 
    
phaselines = 
      Table[phaseline[t], {t, tlow, thigh, \(thigh - tlow\)\/8}]; 


    Show[setup, phaselines]; \)

See the phase lines without the flow plot:

Show[solutionplots, phaselines, Axes -> True, 
      AxesLabel -> {"t", "y"}, PlotRange -> {ylow, yhigh}, 
      AspectRatio -> (1/GoldenRatio), 
      DisplayFunction -> $DisplayFunction];

This gives some (but not all) of the info in the flowplot:

Show[setup];

What a great tool this computer is.

B.3)  Autonomous diffeqs with parameters. Bifurcations and bifurcation points

B.3.a.i)

You are the manager of the Red Oak Catfish Farm in Hartland, Wisconsin. Today your attention is focused on a lake where the logistic model
       y'[t]=a y[t] (1 - y[t]/b) with 0 < a and 0 < b
is used to estimate the fish population t weeks after the lake is stocked. 
Incorporate a constant weekly harvest rate r. The model
       y'[t]=a y[t] (1 - y[t]/b) - r,
is reasonable.
The specifics for this lake are in the code below:

\(Clear[diffeq, y, t, f, r, starter]; 
a\  = \ 0.34; 
b\  = \ 7.1; 

    f[t_, y_]\  = \ a\ y\ \((1 - y\/b)\) - r; 


    setup\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})
          \) // ColumnForm\)

The secrets of this model are contained in the function f[t, y]:

f[t, y]

For each constant harvest rate r, this diffeq is autonomous.
Reason: f[t, y] has no dependence on t.

Take a look at this labeled contour plot of f[t, y] as a function of y and the harvest rate r:

{ylow, yhigh} = {0, 8}; \({rlow, rhigh} = {0, 0.7}; 


    fcontourplot = 
      ContourPlot[f[t, y], {r, rlow, rhigh}, {y, ylow, yhigh}, 
        ContourSmoothing -> Automatic, PlotPoints -> 50, 
        Contours -> {0}, Axes -> True, 
        AxesLabel -> {"r", "y"}, 
        Epilog -> {{Blue, Text["f[t,y] > 0", {0.2, 3}]}, {Orange, 
              Text["f[t,y] < 0", {0.5, 7.5}]}}]; \)

Because 
         y'[t]  =  f[t, y[t]], 
you can relabel this plot:

phaseplot = 
      Show[fcontourplot, 
        Epilog -> {{Blue, Text["y' > 0", {0.2, 3}]}, {Orange, 
              Text["y' < 0", {0.5, 7.5}]}}];

What does this plot tell you ?

Take another look:

Show[phaseplot];

At the surface level, it tells you:
1)  If you have a fish population y and go with harvest rate r so that {r, y} plots out in the black region, then the fish population y is going down. 
2)  If you have a fish population y and go with harvest rate r so that {r, y} plots out in the white region, then the fish population y is going up.

For example, if the harvest rate is r = 0.4 and the fish population is y = 7.0, then the point {r, y} = {0.4, 7.0} plots out as follows:

{r, y} = {0.4, 7.0}Show[phaseplot, Graphics[{Red, PointSize[0.04], Point[{r, y}]}], 
      PlotLabel -> "y going down"];

This signals when you go with harvest rate r = 0.4 and the fish population is y = 7.0, then the fish population is going down. 

See what happens when you go with harvest rate r = 0.4 and the fish population is y = 5.0:

{r, y} = {0.4, 5.0}Show[phaseplot, Graphics[{Red, PointSize[0.04], Point[{r, y}]}], 
      PlotLabel -> "y going up"];

This signals when you go with harvest rate r = 0.4 and the fish population is y = 5.0, then the fish population is going up.

B.3.a.ii)

Take another look at the phase plot:

Show[phaseplot];

Use this plot to say everything you can about how the population y behaves when you go with a constant harvest rate r = 0.4.

Answer:

Throw in a phase line corresponding to r = 0.4:

\(Clear[phaseline, y]; 
r = 0.4; 

    \(phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Red, r]; \)\)newphaseplot = 
      Show[phaseplot, phaseline, 
        PlotRange -> {{rlow, rhigh}, {ylow, yhigh}}];

This tells you that when you maintain a constant harvest rate 
       r = 0.4, 
then: 
→ If you start with y[0] above about 6 (thousand), then the fish population will decrease and level off at about 6 (thousand).
→ If you start with y[0] below about 6 (thousand),  but above about 1.5 , then the fish population will increase and level off at about 6 (thousand).
→ If you start with y[0] below about 1.5 (thousand)  then the fish population will decrease to extinction.
See all this happen:

r = 0.4; endtime = 52; \(Clear[y1, y2, y3, y4, y, t]; 
    
{starter1, starter2, starter3, starter4} = {1.487, 2.0, 4.2, 8.0}; 

    \(y1[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; \)\)y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; y3[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}][[1
          ]]; y4[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter4, 
            y[t], {t, 0, endtime}][[1
          ]]; \(solutionplots = 
      Plot[{y1[t], y2[t], y3[t], y4[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 


    outcomes = 
      Show[solutionplots, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)

Compare:

Show[GraphicsArray[{outcomes, newphaseplot}]];

The new phase plot tells you all you need to know.

B.3.a.iii)

Take another look at the phase plot:

Show[phaseplot];

Use this plot to say everything you can about how the population y behaves when you go with a constant harvest rate r = 0.65.

Answer:

Throw in a phase line corresponding to r = 0.65:

\(Clear[y]; 
r = 0.65; 

    \(phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Red, r]; \)\)newphaseplot = 
      Show[phaseplot, phaseline, 
        PlotRange -> {{rlow, rhigh}, {ylow, yhigh}}];

This tells you that when you maintain a constant harvest rate of
       r = 0.65, 
then the fish population eventually dies out.
See it happen:

r = 0.65; endtime = 52; \(Clear[y1, y2, y3, y4, y, t]; 
    
{starter1, starter2, starter3, starter4} = {1.487, 2.0, 4.2, 8.0}; 

    \(y1[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; \)\)y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; y3[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y4[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter4, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t], y3[t], y4[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 

    \(outcomes = 
      Show[solutionplots, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)\)

Disaster.
Compare:

Show[GraphicsArray[{outcomes, newphaseplot}]];

Again, the new phase plot tells you all you need to know.

B.3.a.iv)

Use the phase plot to estimate the largest sustainable harvest rate r. Estimate how large the starting value y[0] will have have to be to be able to maintain this largest sustainable harvest rate.

Take a peek:

Show[phaseplot];

The largest sustainable harvest rate r is about r ≈ 0.6:

r = 0.6; phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Red, r]; newphaseplot = 
      Show[phaseplot, phaseline, 
        PlotRange -> {{rlow, rhigh}, {ylow, yhigh}}];

And to sustain this harvest rate you'll have to start with y[0] more than about 3.3 (thousand fish). Check it out:

r = 0.6; endtime = 52; \(Clear[y1, y2, y3, y4, y, t]; 
    
{starter1, starter2, starter3, starter4} = {2.9, 3.4, 4.2, 8.0}; 

    \(y1[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; \)\)y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; y3[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y4[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter4, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t], y3[t], y4[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 


    outcomes = 
      Show[solutionplots, Axes -> True, 
        AxesLabel -> {"t", "y"}, 
        PlotRange -> {ylow, yhigh}, AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)

Compare:

Show[GraphicsArray[{outcomes, newphaseplot}]];

B.3.a.v) Estimating the bifurcation point

Take yet another look at the phase plot this time shown with a line indicating the largest sustainable harvest rate estimate:

r = 0.6; bifurcationplot = 
      Show[phaseplot, Graphics[{Red, Line[{{r, ylow}, {r, yhigh}}]}]];

Lots of folks call this a bifurcation plot.  Neat word.
Why do they call it by that funny name?
Estimate the birfurcation point.

A birfurcation is any situation in which you see something forking into two branches. 
And you are seeing it here:

Show[bifurcationplot];

Now look at the diffeq again;

\(Clear[r]; 
diffeq\)

If you make r bigger than the r indicated in the plot, then ALL solutions of the diffeq decay to 0.
If you make r lower than the r indicated in the plot, then some of the solutions decay to 0, but others do not.

The cutoff point between the differences in solution behavior is the r indicated in the plot.

That's why this plot is called a bifurcation plot. And that's why the r indicated in the plot is called the bifurcation point.

B.4)  Sensitive dependence on starter data

B.4.a.i) Flow plots

Here's a simple differential equation:

\(Clear[diffeq, y, t, f, starter]; 
b\  = \ 10; 

    f[t_, y_]\  = \ y\  - \ 0.8\ t; 


    setup\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})
          \) // ColumnForm\)

And a flow plot together with the plots of two solutions::

endtime = 6; scalefactor = 0.4; {ylow, yhigh} = {-4, 10}; \({tlow, thigh} = {0, 10}; 


    flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.7], {t, tlow, thigh, \(thigh - tlow\)\/10}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/8}]; 

Clear[y1, y2, y3, t]; 
    
{starter1, starter2} = {0.79, 0.81}; 
\(endtime = 7; \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y2[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, 
        DisplayFunction -> Identity]; 


    Show[solutionplots, flowplot, DisplayFunction -> $DisplayFunction]; 
    \)

What do folks mean when they talk about "sensitive dependence on starter data?"

You're looking at it.  
Take another look:

Show[solutionplots, flowplot, 
      DisplayFunction -> $DisplayFunction];

What you see are two solutions of:

diffeq

The starting points of the two solutions are almost the same; yet the two solutions are totally different in character. 
You gotta agree that small changes in the starting conditions can result in big changes in solutions. This is exactly what folks mean when they talk about senistive dependence on starting data.

B.4.a.ii)

Do all differential equations exhibit sensitive dependence on starter data?

Not all of them.
Here's one that doesn't exhibit sensitive dependence on starter data:

\(Clear[f, t, y]; 
b\  = \ 10; 
f[t_, y_]\  = \ \ \(-y\)\ ; 


    setup\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})
          \) // ColumnForm\)

endtime = 6; scalefactor = 0.4; {ylow, yhigh} = {-4, 4}; \({tlow, thigh} = {0, 8}; 


    flowplot = 
      Table[Arrow[scalefactor\ {1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue], {t, tlow, thigh, \(thigh - tlow\)\/8}, {
          y, ylow, yhigh, \(yhigh - ylow\)\/6}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"t", "y"}, 
      AspectRatio -> 1\/GoldenRatio, 
      PlotLabel -> setup[[1, 1]]]; 
    \)

The flow plot tells you to expect all solutions to the sucked onto the t-axis.
Here it is for three solutions:

\(Clear[y1, y2, y3, t]; 
    
{starter1, starter2, starter3} = {\(-2.0\), 2.0, 4.0}; 

    endtime = thigh; 


    y1[t_] = y[t] /. 
        \(NDSolve[diffeq /. starter -> starter1, y[t], {t, 0, endtime}]
            \)[[1]]; 

    \(y2[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; \)\)\(y3[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, 
        DisplayFunction -> Identity]; 


    Show[solutionplots, flowplot, 
      PlotLabel -> setup[[1, 1]], 
      DisplayFunction -> $DisplayFunction]; \)

If two solutions of this differential equation have nearby starting values on y[0], then the corresponding solutions will remain close to each other forever.
This differential equation is nearly insensitive to starter data.

B.4.b) Using the phase line to detect extreme sensitivity to errors in starter data

Here's an autonomous diffeq:

\(Clear[diffeq, y, t, f, starter]; 
b\  = \ 6; 
c\  = \ 2; 

    f[t_, y_]\  = \ 0.31\ y\ \((1 - y\/b)\)\ \((1 - y\/c)\); 


    \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\) // 
      ColumnForm\)

Here is a phaseline for this diffeq:

{ylow, yhigh} = {0, 10}; \({tlow, thigh} = {0, 15}; 


    phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Blue, \(-1\)]; 


    setup = Show[phaseline, AspectRatio -> 1\/GoldenRatio, 
        Axes -> True, 
        PlotRange -> {{tlow - 2, thigh}, {ylow, yhigh}}, 
        AxesLabel -> {"t", "y[t]"}]; \)

Use the phase line to determine starting values on y[0] at which this diffeq is extremely sensitive to error in y[0].

Take another look:

Show[setup];

If y[0] is just below 6, then solutions go down and level off at  y = 2.
If y[0] is just above 6, then solutions go up.
So you have sensitive dependence on starter data for y[0] near 6.
Check it out:

endtime = thigh; \({starter1, starter2} = {5.95, 6.05}; 


    y1[t_] = y[t] /. 
        \(NDSolve[diffeq /. starter -> starter1, y[t], {t, 0, endtime}]
            \)[[1]]; 

    y2[t_] = y[t] /. 
        \(NDSolve[diffeq /. starter -> starter2, y[t], {t, 0, endtime}]
            \)[[1]]; 


    solutionplots = 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, 
        DisplayFunction -> Identity]; 


    Show[setup, solutionplots, DisplayFunction -> $DisplayFunction]; 
    \)

Now that's sensitive!

B.4.c) Using bifurcation plots to detect extreme sensitivity to errors in starter data

You are the manager of the Red Oak Catfish Farm in Hartland, Wisconsin. Today your attention is focused on a lake where the logistic model
    y'[t]=a y[t] (1 - y[t]/b) with 0 < a and 0 < b
is used to estimate the fish population t weeks after the lake is stocked. 
Incorporate a constant weekly harvest rate r. The model
    y'[t]=a y[t] (1 - y[t]/b) - r,
is reasonable.
The specifics for this lake are in the code below:

\(Clear[diffeq, y, t, f, r, starter]; 
a\  = \ 0.34; 
b\  = \ 7.1; 

    f[t_, y_]\  = a\ y\ \((1 - y\/b)\) - r; 


    setup\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})
          \) // ColumnForm\)

The secrets of this model are contained in the function f[t, y]:

f[t, y]

For each constant harvest rate r, this diffeq is autonomous.
Reason: f[t, y] has no dependence on t.

Take a look at this labeled contour plot of y'=f[t, y] as a function of y and the harvest rate r:

{ylow, yhigh} = {0, 8}; \({rlow, rhigh} = {0, 0.7}; 


    phaseplot = 
      ContourPlot[f[t, y], {r, rlow, rhigh}, {y, ylow, yhigh}, 
        ContourSmoothing -> Automatic, PlotPoints -> 50, 
        Contours -> {0}, Axes -> True, 
        AxesLabel -> {"r", "y"}, 
        Epilog -> {{Blue, Text["y' > 0", {0.2, 3}]}, {Orange, 
              Text["y' < 0", {0.5, 7.5}]}}]; \)

What does this plot tell you about sensitive dependence on starter conditions for this diffeq?

Take another look:

Show[phaseplot];

This also tells you that when you go with a harvest rate 
         r less than about 0.6, 
then you should expect sensitive dependence to any starter data on y[0] that is near the bottom half of the curve. 
If the starter data on y[0] is a little bit above the bottom half of the curve, then the corresponding solution will increase but if the starter data on y[0] is just a little bit below the curve, then the corresponding solution will decrease.

r = Random[Real, {rlow, 0.58}]; \(phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Red, r]; 


    Show[phaseplot, phaseline, 
      PlotRange -> {{rlow, rhigh}, {ylow, yhigh}}]; \)

r = Random[Real, {rlow, 0.58}]; \(phaseline = PhaseLine[f[t, y], {y, ylow, yhigh}, Red, r]; 


    Show[phaseplot, phaseline, 
      PlotRange -> {{rlow, rhigh}, {ylow, yhigh}}]; \)

Tremendous sensitivity for {r, y[0]} near the bottom half of the curve. 
No serious sensitivity for {r, y[0]} along the top half of the curve.
And no serious sensitivity elsewhere.

B.5)  Formulas for solutions: Nice work when you can get it

Rensselaer Polytechnic Institute mathematics professor William Boyce and  co-author of the respected standard DiffEq textbook "Elementary Differential Equations and Boundary Value Problems (Wiley, 1986) says this:

Stanford University electrical engineering professor Steve Boyd puts it this way:

Stanford University electrical engineering professor and National Academy of Engineering member, Tony Siegman says this:

B.5.a.i)

If you have a formula for f[t, y], then you can always plot the flow of solutions of
      y'[t]=f[t, y[t]].
 For example:

\(Clear[f, t, y]; 
f[t_, y_]\  = \ \ 0.3\ y\ \((y\  - \ \ 0.4\ t)\); 

    diffeq\  = \ \((\(y'\)[t]\  == \ f[t, y[t]])\)\)

And a flow plot:

scalefactor = 0.5; {ylow, yhigh} = {0, 6}; \({tlow, thigh} = {0, 10}; 


    flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> Normalize, 
          HeadSize -> 0.5], {t, tlow, thigh, \(thigh - tlow\)\/8}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/6}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"t", "y"}, 
      AspectRatio -> 1\/2]; \)

And you can plot a sample solution:

\(Clear[y, ysol, t]; 
endtime = 12; 
thisstarter = 0.9; 


    sol = NDSolve[{diffeq, y[0] == thisstarter}, y[t], {t, 0, endtime}]; 

    ysol[t_] = y[t] /. sol[[1]]; 


    solutionplot = 
      Plot[ysol[t], {t, 0, endtime}, PlotRange -> All, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, AspectRatio -> 1\/2, 
        DisplayFunction -> Identity]; 


    Show[solutionplot, flowplot, DisplayFunction -> $DisplayFunction]; 
    \)

Say, it looks like this particular solution goes to  0 as t -> ∞.
How can you tell for sure?

You ask Mathematica for a formula for the solution:

formulasol = DSolve[{diffeq, y[0] == thisstarter}, y[t], t]

Ugh. That's  a mess. 
To get a useful formula, look at the diffeq again:

diffeq

Re-enter it with cleared symbols:

\(Clear[a, b, f, t, y, starter]; 
f[t_, y_]\  = a\ y\ \ \((y - b\ t)\); 
    
diffeq\  = \ \((\(y'\)[t]\  == \ f[t, y[t]])\)\)

Ask Mathematica  for a solution this way:

thisa = 0.3; thisb = 0.4; 
    formulasol = 
      DSolve[{diffeq, y[0] == starter}, y[t], t] /. {a -> thisa, 
          b -> thisb, starter -> thisstarter}

Fish out the formula:

\(Clear[yformula]; 

    yformula[t_] = 
      Chop[y[t] /. formulasol[[1]]]\)

As t -> ∞, Erf[t] -> 1:

Plot[Erf[t], {t, 0, 8}, PlotStyle -> {{Thickness[0.01], Blue}}, 
      PlotRange -> All, AspectRatio -> (1/GoldenRatio), 
      AxesLabel -> {"t", "Erf[t]"}];

So as t -> ∞, the dominant term is the (e^(-0.06 (t^2))) term which goes to 0 as  t -> ∞. This forces the solution to go to 0 as t -> ∞:

Limit[yformula[t], t -> ∞]

Now you know for sure that the solution plotted above is 0 in the global scale.

B.5.a.ii)

If you have a formula for f[t, y], then you can always plot the flow of solutions of
      y'[t]  =  f[t, y[t]].
 For example:

\(Clear[f, t, y]; 
f[t_, y_]\  = \ \ 0.3\ y\ Sin[0.4\ \((y\  + \ t)\)]; 
    
diffeq\  = \ \((\(y'\)[t]\  == \ f[t, y[t]])\)\)

And a flow plot:

scalefactor = 0.8; {ylow, yhigh} = {0, 20}; \({tlow, thigh} = {0, 20}; 


    flowplot = 
      Table[Arrow[scalefactor\ {1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> 1, 
          HeadSize -> 0.6], {t, tlow, thigh, \(thigh - tlow\)\/15}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/15}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"t", "y"}, 
      AspectRatio -> 1\/GoldenRatio]; \)

And you can plot a sample solution:

\(Clear[y, ysol, t]; 
endtime = thigh; 
starter = 16; 


    sol = NDSolve[{diffeq, y[0] == starter}, y[t], {t, 0, endtime}]; 

    ysol[t_] = y[t] /. sol[[1]]; 


    solutionplot = 
      Plot[ysol[t], {t, 0, endtime}, PlotRange -> All, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        AxesLabel -> {"t", "y[t]"}, AspectRatio -> 1\/2, 
        DisplayFunction -> Identity]; 


    Show[solutionplot, flowplot, DisplayFunction -> $DisplayFunction]; 
    \)

It looks like this particular solution gets sucked onto a straight line.
How can you tell for sure?

You could ask Mathematica for a formula for the solution, but Mathematica will fail. This is not because of any deficiency of Mathematica. 

No machine, no person, no anything has ever come up with a formula for the solution of this differential equation. 
In this case, you'll have to depend on the visual evidence because the confirming formula is simply not available.

B.5.a.iii) Formulas are luxuries

What's the moral?

Professor Peter Lax, former President of the American Mathematical Society, once said that the class of differential equations whose solutions can be described by exact formulas is "pitifully small."

Students in differential equations courses of the past have been forced to concentrate on this small class to the exclusion of others. You have bigger fish to fry.

Think of it this way:
If you can get a neat simple formula for a solution of a given differential equation, then you are very lucky. 

According to Lax, the main thing to remember is "the general idea that every differential equation has a solution and that this solution is uniquely determined by initial data. . . That today we can use computers to explore the solutions of [differential] equations is truly revolutionary; we are only beginning to glimpse the consequences."

You are in on the ground floor of this revolution.