DE.05 First Order Differential Equations
TUTORIALS

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];

T.1)  Setting harvest rates to control the catfish population

T.1.a.i)

You are the manager of the Black Oak Catfish Farm in Moberly, Missouri. One day a team from the University of Missouri College of Agriculture stops by with some advice based on new research. First, they tell you that pushing the population too high invites disease and death of fish.  They go on to say that you should try to set your harvest rate so that the long term fish population is 0.7 times the carrying capacity of the lake.

You decide to apply this to a certain lake modeled by 
          y'[t]=a y[t] (1 - y[t]/b) - r 
where y[t] measures the population of the lake (in thousands), a = 0.12, b = 9 (thousand) and r is your weekly harvest rate (in thousands).
Come up with an estimate of the r that does this.

Knowing that a bifurcation diagram reveals a heckuva lot, you whip out your PowerPC laptop, ignite Mathematica and set up your model:

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

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


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

You know that the secrets of the model are contained in:

f[t, y]

You make this labeled contour plot of f[t, y] as a function of y and the weekly harvest rate r:

{ylow, yhigh} = {0, b}; \({rlow, rhigh} = {0, 0.35}; 


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

The carrying capacity of the lake is b and you want to control the long term population to be 0.7 b,  so you plot:

ygoal = 0.7 bgoalplot = 
      Show[phaseplot, 
        Graphics[{Green, Thickness[0.01], 
            Line[{{rlow, ygoal}, {rhigh, ygoal}}]}], 
        PlotLabel -> "Goal plot"];

You want 
         y'[t]  =  0 when y[t] = ygoal 
(so that the population doesn't change when it reaches ygoal), so you want to set r so that {r, ygoal} is the point at which the goal line goes from the light region (where y'[t]>0) to the black region (where y'[t]<0). 
Because 
         y'[t]=f[t, y[t]], 
the r you want is the r that makes
         f[t,  ygoal] = 0:

\(Clear[r]; 
controleqn = 0 == f[t, y[t]] /. y[t] -> ygoal; 

    controlsol = Solve[controleqn, r]; 


    rcontroller = 
      r /. controlsol[[1]]\)

Or:

FindRoot[controleqn, {r, 0.2}]

This is the harvest rate the model tells you to use.
Take a look:

\(controlline = 
      Graphics[{Green, Thickness[0.01], 
          Line[{{rcontroller, ylow}, {rcontroller, yhigh}}]}]; 


    Show[goalplot, controlline, 
      Graphics[{Green, PointSize[0.04], Point[{rcontroller, ygoal}]}]]; \)

Try out this harvest rate r for several starting values on y[0]:

r = rcontroller; \(Clear[y1, y2, y3, t]; 
{starter1, starter2, starter3} = {3, 5, 8}; 

    \(endtime = 150; \)\)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
          ]]; 


    Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Thickness[0.015], Red}}, 
      AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
      AxesLabel -> {"t", "y[t]"}, 
      Epilog -> {Green, Thickness[0.01], 
          Line[{{0, ygoal}, {endtime, ygoal}}]}]; \)

Yes! 
The long term population seems to be controlled just as it was supposed to be.

T.1.a.ii)

Stay with the model in part i) above, and make sure all parts of part i) are live on your machine and look at this:

r = rcontroller; \(Clear[y1, y2, y3, t]; 
    
{starter1, starter2, starter3} = {2.67, 3.0, 5.0}; 
\(endtime = 150; 
    \)\)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
          ]]; 


    Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Thickness[0.015], Red}}, 
      AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
      AxesLabel -> {"t", "y[t]"}, 
      Epilog -> {Green, Thickness[0.01], 
          Line[{{0, ygoal}, {endtime, ygoal}}]}]; \)

Apparently setting the harvest rate
        r = rcontroller 
works for some starting values on y[0] and not for others. 

How big does the starting value on y[0] have to be to guarantee that setting
        r = rcontroller
does the job it was selected to do?

Look at:

startplot = 
      Show[phaseplot, goalplot, 
        Graphics[{Yellow, Thickness[0.01], 
            Line[{{rcontroller, ylow}, {rcontroller, yhigh}}]}]];

If you go with a starter value on y[0] so that {rcontroller, y[0]} lands on the plotted line in the lower part of the black region, then you are in trouble because you are starting with y'[0]<0 and things will only get worse as time goes on. 
You can estimate from the plot that 
         r = rcontroller 
will not work for starting values on 
         y[0] < 2.5 
and that it will work fine for starting values on 
         y[0] > 3.0:

r = rcontroller; \(Clear[y1, y2, t]; 
{starter1, starter2} = {2.5, 3.0}; 
endtime = 100; 
    

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
          ]]; 


    twosolutions = 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, 
        Epilog -> {Green, Thickness[0.01], 
            Line[{{0, ygoal}, {endtime, ygoal}}]}]; \)\)

You can refine your estimate of the critical starting value by plotting:

\(Clear[y1, y2, t]; 
{starter1, starter2} = {2.6, 2.7}; 
endtime = 100; 
    

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]]; 


    twosolutions = 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, 
        Epilog -> {Green, Thickness[0.01], 
            Line[{{0, ygoal}, {endtime, ygoal}}]}]; \)

The flow plot confirms that when you go with y[0] > 2.7, then the control will work, but when you go with y[0] < 2.6, then the control will not work. With further experimentation, you can refine this estimate.

T.1.a.iii)

The University of Missouri team comes back to your fish farm and you excitedly tell them what you've learned.  They are happy to see your analysis.  They go on to say, "The model
        y'[t]=a y[t] (1 - y[t]/b) - r
treats every day the same way.  Bright summer days make for faster growth and higher lake capacity than the dark days of winter."
In fact, the proportion of the day getting daylight is given by this function (with t measuring weeks since January 1):"

\(Clear[daylightproportion, t]; 

    daylightproportion[t_] = 
      1\/24\ \((12.00 + 2.4\ Sin[0.0172\ \((\(365\ t\)\/52 - 80)\)])\)\)

You plot this function:

Plot[daylightproportion[t], {t, 0, 52}, 
      PlotStyle -> {{Thickness[0.01], Red}}, AspectRatio -> (1/4)];

You say: "This looks good because it's high in the summer and low in the winter. Now let me check it by calculating 
        (1/52) 
        (∫_0^52) daylightproportion[t] dt:"

(1/52) NIntegrate[daylightproportion[t], {t, 0, 52}]

And you react: "I like it because now I know for sure that the average value of daylightproportion[t] is very close to 0.5, as it should be."
The University of Missouri teams says that the growth rate a varies with sunlight as does the lake capacity b. And they say that you should refine the model by replacing the constant growth factor a by a function whose average value is a,  but is proportional to the hours of daylight on day t and whose average value is 1 like this:

\(Clear[bettera]; 
bettera[t_] = a\ 2\ daylightproportion[t]\)

And you should replace the constant lake capacity b by b times the same function:

\(Clear[betterb]; 
betterb[t_] = b\ 2\ daylightproportion[t]\)

You say: "Let's see what this refined model predicts for two years with a weekly harvest rate r  =  0.2 (thousand) fish but starting out with 4000 fish."

\(r\  = \ 0.2; 
Clear[diffeq, y, t, f, starter]; 


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

    model\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})\); 
    
diffeq // ColumnForm
\)\(Clear[y1, y2, t]; 
{starter1} = {4}; 
\(endtime = 2\ 52; \)\)ygoal = 6.3; \(y1[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


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


    goalplot = 
      Graphics[{Green, Thickness[0.01], 
          Line[{{0, ygoal}, {endtime, ygoal}}]}]; 
\(scalefactor = 4.5; 
    \)\)\({ylow, yhigh} = {0, 9}; 


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


    Show[flowplot, goalplot, solutionplot, 
      AspectRatio -> 1\/GoldenRatio, Axes -> True, 
      AxesLabel -> {"t", "y[t]"}, PlotRange -> All, 
      DisplayFunction -> $DisplayFunction]; \)

You say: "I like this.  Now let me try to use a variable harvest rate r[t] to control the fish population so that, in the long term, y[t] settles in on the goal of 6.3 (thousand) fish."
Do it. Estimate the size of the first year's harvest.

You do the same thing you did above but incorporate a variable weekly harvest rate r[t]:

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


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


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

You want to maintain the fish population at 6.3 (thousand); so
           y'[t]  =  f[t, y[t]]
should be 0 when y[t] =  6.3:

controleqn = 0 == f[t, y[t]] /. y[t] -> 6.3

Your variable weekly harvest rate r[t] is:

sol = Solve[controleqn, r[t]]; r[t_] = r[t] /. sol[[1]]

Check it out:

\(Clear[y1, y2, y3, t]; 
    
{starter1, starter2, starter3} = {3.0, 5.0, 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
          ]]; 


    solutionplot = 
      Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        DisplayFunction -> Identity]; 
\(scalefactor = 3.5; \)\)\({ylow, yhigh} = {0, 9}; 


    flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> 6, 
          HeadSize -> 0.6], {t, 0, endtime, endtime\/10}, {y, ylow, 
          yhigh, \(yhigh - ylow\)\/8}]; 


    Show[flowplot, goalplot, solutionplot, 
      AspectRatio -> 1\/GoldenRatio, Axes -> True, 
      AxesLabel -> {"t", "y[t]"}, PlotRange -> All, 
      DisplayFunction -> $DisplayFunction]; \)

Looks good and feels good.
In all these cases, the first year's harvest brings in this many fish:

1000 NIntegrate[r[t], {t, 0, 52}]

Not bad.

T.1.a.iv)

Keep everything (except the starting values on y[0]) the same as in the previous part and estimate the minimum starter value on y[0] that will guarantee that the fish population will not die out.

Take a look at the flow plot:

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

It looks like the minimum starter value of y[0] that will guarantee that the fish population will not die out is around y[0] ≈ 3.0.
Check it out:

\(Clear[y1, y2, y3, t]; 
    
{starter1, starter2, starter3} = {2.5, 3.0, 3.5}; 

    \(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
          ]]; 


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


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

The three starter values in the plot above were {2.5, 3.0, 3.5}, and the starter value of 2.5 did not lead to a sustainable harvest. This leads to the estimate y[0]  ≈  3.0 for the smallest starter value for a sustainable harvest.
You can refine this estimate :

\(Clear[y1, y2, y3, t]; 
    
{starter1, starter2, starter3} = {2.7, 2.9, 3.1}\)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
          ]]; 


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


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

The three starter values for this plot were: {2.7, 2.9  .3  .1}. 
A good estimate of the minimum starter value on y[0] that will guarantee that the fish population will not die out is 
          y[0] ≈ 2.7 
but when you start with this value, you are playing with fire.

T.2)  Hand symbol manipulation: Separating the variables

T.2.a)

All oldtime DiffEq folks have heard of the method of separation of variables. 
Just what is this famous procedure?

It's a way of using symbolic manipulation to go after an exact formula for the solution of certain differential equations. It is normally used only in simple, special cases.
Here it is in action on the following diffeq:

\(Clear[diffeq, y, t, f, r, starter]; 
f[t_, y_]\  = \(-2.1\)\ y\^2; 


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

Before Mathematica will help you, you have to re-enter the diffeq with cleared constants this way:

\(thisa\  = \ \(-2.1\); 
thisstarter\  = \ 3.1; 

    Clear[a, y, t, diffeq, y, t, f, r, starter]; 

    f[t_, y_]\  = \ \ a\ y^2\ t; 


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

Separating the variables involves putting all the y[t] and y'[t] terms on the left and leaving the t terms on the right like this:

newleft=((y'[t])/(y[t]^2));newright = a t; samediffeq = newleft == newright

Now integrate both sides from 0 to t:

integratedeqn = 
      (∫_0^t) newleft \[DifferentialD]t == 
        (∫_0^t) newright \[DifferentialD]t

Replace y[0] by  starter:

neweqn = integratedeqn /. y[0] -> starter

Now solve for y[t]

\(Clear[symbolformulay]; 
ysol = Solve[neweqn, y[t]]; 

    symbolformulay[t_] = 
      y[t] /. ysol[[1]]\)

Now replace a with the given a and starter with the given starter:

\(\(Clear[formulay]; 

    formulay[t_]\  = \ 
      symbolformulay[t] /. {a -> thisa, starter -> thisstarter}\ \)\)

Done.

T.2.b)

Here's another diffeq:

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

    f[t_, y_]\  = t\ \((1.2 + 2.4\ y\^2)\); 


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

Use separation of variables to try to come up with a formula for the solution.

Take another look:

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

    f[t_, y_]\  = t\ \((1.2 + 2.4\ y\^2)\); 


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

Before Mathematica will help you, you have to re-enter the diffeq with cleared constants this way:

\(thisa\ \  = \ 1.2; 
thisb\  = \ 2.4; 

    Clear[a, y, t, diffeq, y, t, f, r, starter]; 

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


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

Separating the variables involves putting all the y[t] and y'[t] terms on the left and leaving the t terms on the right like this:

newleft=((y'[t])/(a + b y[t]^2));newright = t; samediffeq = newleft == newright

Now integrate both sides from 0 to t:

integratedeqn = 
      (∫_0^t) newleft \[DifferentialD]t == 
        (∫_0^t) newright \[DifferentialD]t

Replace y[0] by  starter:

neweqn = integratedeqn /. y[0] -> starter

Now solve for y[t]

\(Clear[symbolformulay]; 
ysol = Solve[neweqn, y[t]]; 

    symbolformulay[t_] = 
      y[t] /. ysol[[1]]\)

Now replace a with the given a and b with the given b:

\(\(Clear[formulay]; 

    formulay[t_]\  = \ symbolformulay[t] /. {a -> thisa, b -> thisb}\ \)\)

Done.

T.2.c) The logistic diffeq

Here's the logistic diffeq:
          y'[t]=r y[t] (1 - y[t]/b)
          with y[0] = starter.

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

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


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

Here's Mathematica  cranking out a formula for the solution:

DSolve[diffeq, y[t], t]

Use separation of variables to explain where this formula comes from.

Look at the diffeq again:

diffeq

Separate the variables:

newleft=((y'[t])/(y[t] (1 - y[t]/b)));newright = r; samediffeq = newleft == newright

Integrate both sides from 0 to t:

integratedeqn = 
      (∫_0^t) newleft \[DifferentialD]t == 
        (∫_0^t) newright \[DifferentialD]t

Replace y[0] with the (cleared) starter value on y[0]:

neweqn = integratedeqn /. y[0] -> starter

And solve for y[t]:

\(Clear[formulay]; 
ysol = Solve[neweqn, y[t]]\)

There it is.

T.2.d) Sometimes separation of variables is useless

Here are two new differential equations:

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

    f[t_, y_]\  = \ \ \(-2.1\)\ \ Sin[\ y\  + \ t]\  + \ y; 


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

\(Clear[diffeq, y, t, f, r, starter]; 
f[t_, y_]\  = y\^2 + t; 


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

Why is separation of variables useless for these two differential equations?

Look at the first one:

diffeq1

Try as you may, you can't separate the y[t]'s from the t's.

Look at the second:

diffeq2

This time you can separate the variables:

newleft=y'[t]-y[t]^2;newright = t; samediffeq = newleft == newright

But when you integrate both sides from 0 to t, you get:

integratedeqn = 
      (∫_0^t) newleft \[DifferentialD]t == 
        (∫_0^t) newright \[DifferentialD]t

The integral on the left could not be done.
Mathematica  has its own obscure ways of solving this:

DSolve[diffeq2, y[t], t]

But the formula is overwhelming.
To get control of this one, go to flow plots and NDSolve.

T.3)  The danger zone: Where the unexpected sometimes happens and separation of variables produces bogus formulas

T.3.a.i)

Here is an apparently innocent-looking differential equation:

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

    f[t_, y_]\  = \ \ 0.26\ t\ \((\@\(1 - y\^2\))\); 


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

You can try to generate your own formula for the solution of this differential equation by separating the variables:

Separate and integrate:

starter  =  0.2;
a  =  0.26;
Clear[symbolyformula];
newleft=((y'[t])/((Sqrt[1 - y[t]^2])));newright = a t; \(neweqn = 
      \((∫\_0\%t newleft \[DifferentialD]t /. y[0] -> starter)
          \) == ∫\_0\%t newright \[DifferentialD]t; 

    ysol = Solve[neweqn, y[t]]; 


    symbolformulay[t_] = \ 
      N[ExpandAll[y[t] /. ysol[[1]]], 6]
        \)

This is the same as the formula delivered by DSolve:

dsol=DSolve[{y'[t]==f[t, y[t]],y[0] == 0.2},y[t],t];N[ExpandAll[y[t] /. dsol[[1]]], 6]

Here is the way this formula plots out:

\(endtime\  = \ 5; 

    formulaplot\  = \ 
      Plot[symbolformulay[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AxesLabel -> {"[t", "y[t]"}, AspectRatio -> 1/3]; \)

See this plot along with the flow plot for the given diffeq:

{ylow, yhigh} = {-0.2, 1.0}; \({tlow, thigh} = {0, endtime}; 


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


    badplot = 
      Show[formulaplot, solutionflowplot, Axes -> True, 
        AxesLabel -> {"t", "y[t]"}, 
        AspectRatio -> 1\/GoldenRatio, PlotLabel -> "Bad"]; 
    \)

The solution plot goes with the flow on the left but not on the right. 
On the right, this plot is dead wrong.
Fix it.

This plot is also out of sync with the differential equation:

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

    f[t_, y_]\  = \ \ 0.26\ t\ \((\@\(1 - y\^2\))\); 


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

Because 
        0.26 t (Sqrt[1 - y[t]^( 2)]) ≥ 0 for t ≥ 0 
and because
         y'[t]=0.26 t (Sqrt[1 - y[t]^( 2)])for t ≥ 0, 
no solution of this differential equation can ever go down for t ≥ 0.  
This fact is key to understanding what the correct formula must be.

Take another look:

hitsol = FindRoot[symbolformulay[t] == 1, {t, 3}]

marker = Graphics[{PointSize[0.07], Yellow, Point[{3.24, 1}]}]; 
    Show[badplot, marker];

And look at the differential equation again:

diffeq

When y[t] < 1, it's automatic that y[t] goes up because 
        y'[t]=(Sqrt[1 - y[t]^( 2)]) > 0.
But it's impossible for y[t] to be  >  1.
Reason: If y[t] > 1, then (Sqrt[1 - y[t]^( 2)]) is an imaginary number.
The upshot: 
The true solution y[t] steadily increases until y[t] gets to 1 at t = 3.24519. 
For  t > 3.24, 
the true solution y[t] stays equal to 1 because y[t] cannot be bigger than 1 and y[t] cannot go down. 
Thus, for  t > 3.24,  y[t] is trapped equal to 1. 

Watch the true solution go with the flow:

symbolformulay[t]

\(b\  = \ 3.24; 
Clear[truesolution]; 

    truesolution[t_]\  := \ 1 /; t\  > \ b\ ; \ 

    truesolution[t_]\  := \ symbolformulay[t] /; t\  <= \ b

\)\(truesolutionplot\  = \ 
      Plot[truesolution[t], {t, 0, endtime}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle -> {{Thickness[0.015], Red}}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"t", "y[t]"}, DisplayFunction -> Identity]; 

    \ \ \ \ \ \ \ \ \ \ \ \ \ 

    goodplot\  = \ 
      Show[truesolutionplot, solutionflowplot, Axes -> True, 
\ \ \ \ \ 
        AxesLabel -> {"t", "y[t]"}, AspectRatio -> 1\/GoldenRatio, 

        \ \ \ \ \ PlotLabel -> "Good", 
        DisplayFunction -> $DisplayFunction]; \)

Compare:

Show[badplot, marker];

T.3.a.ii) The danger zone and the gradient test for the danger zone

Is there a way of getting an advance alert to sticky situations like this one?

You bet your sweet bippy.
There is a way and it's very simple:

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

    f[t_, y_]\  = \ \ 0.26\ t\ \((\@\(1 - y\^2\))\); 


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

You just look at the gradient of f[t, y]:

\(Clear[gradf]; 

    gradf[t_, y_] = {∂\_t f[t, y], ∂\_y f[t, y]}\)

Notice gradf[t, y] has a singularity (blow up) at y = 1 and y = -1.

gradf[t, 1]

gradf[t, -1]

Yepper. 
This tells you that the possible danger zone for this differential equation consists of all points on the line 
        y[t] = 1 
together with all points on the lines 
        y[t] = -1.
Take a gander:

{ydanger1, ydanger2} = {-1.0, 1.0}; {tlow, thigh} = {0, endtime}; dangerzone = {
        Graphics[{Yellow, Thickness[0.01], 
            Line[{{tlow, ydanger1}, {thigh, ydanger1}}]}], 
        Graphics[{Yellow, Thickness[0.01], 
            Line[{{tlow, ydanger2}, {thigh, ydanger2}}]}]}; dangerlabel = {Graphics[Text["DANGER", {(endtime/2), ydanger1}]], 
        Graphics[Text["DANGER", {(endtime/2), ydanger2}]]}; scalefactor = 0.25; \(solutionflowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.15], {t, tlow, thigh, \(thigh - tlow\)\/8}, {y, 
          ydanger1, ydanger2, \(ydanger2 - ydanger1\)\/8}]; 


    dangerflow = 
      Show[solutionflowplot, dangerzone, dangerlabel, Axes -> True, 
        AxesLabel -> {"t", "y[t]"}, 
        AspectRatio -> 1\/GoldenRatio]; \)

Compare:

Show[dangerflow, truesolutionplot, marker, 
      PlotLabel -> "Good"];

Show[dangerflow, formulaplot, marker, PlotLabel -> "Bad"];

Anytime you see a hand or computer-generated formula plot out and actually touch a danger zone, you should become very alert to the possiblity of having to do a careful analysis.

T.3.a.iii)

Here's a diffeq:

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

    f[t_, y_]\  = \ \(-0.4\)\ \((\@y)\)\ ; 


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

Plot the danger zone for this diffeq.

Calculate the gradient of f[t, y]:

\(Clear[gradf]; 

    gradf[t_, y_] = {∂\_t f[t, y], ∂\_y f[t, y]}\)

Note that gradf[t, y] blows up for y = 0.

gradf[t, 0]

The singularities of gradf[t, y] are all on the line 
           y = 0.
This curve is the danger zone for this diffeq.
Here comes a plot:

endtime = 6; {ylow, yhigh} = {0, 4}; {tlow, thigh} = {0, endtime}; dangerzone = 
      Plot[0, {t, 0, endtime}, PlotStyle -> {{Thickness[0.01], Yellow}}, 
        DisplayFunction -> Identity]; dangerlabel = Graphics[Text["DANGER", {(endtime/4), 0}]]; scalefactor = 0.6; \(flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.3], {t, tlow, thigh, \(thigh - tlow\)\/8}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/8}]; 


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

Anytime you see a solution hit this danger zone, you have to be on guard.

T.3.a.iv)

Here's another diffeq

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

    f[t_, y_]\  = \ 0.4\ Sin[t]\  + \ y\^2\ ; 


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

Plot the danger zone for this diffeq.

Calculate the gradient of f[t, y]:

\(Clear[gradf]; 

    gradf[t_, y_] = {∂\_t f[t, y], ∂\_y f[t, y]}\)

Note that gradf[t, y] never blows up for finite t or y.
This tells you that this diffeq has no danger zone. So there is nothing to plot.
This differential equation can't reach out and surprise you.

T.3.a.v)

Try to explain why the gradient test for the danger zone works

Go with a cleared diffeq:

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

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

Because 
       y'[t]=f[t, y[t]], 
the chain rule from vector calc guarantees that 
       y''[t] = gradf[t,  y[t]]•{1,y'[t]}.

Now here's the key point:
If {t, y[t]} hits the danger zone (i.e. a singularity of gradf[t, y[t]]) at a time t = tdanger, 
then the fact that 
         (y^′)′[tdanger]=gradf[tdanger,  y[tdanger]]•{1,y'[tdanger]}
 opens up the possibility of a singularity of 
         y''[t] at  t = tdanger.
When this happens things can get out of hand! 
And they did in part i) above. In part i), the singularity in  (y^′)′[t] at 
         tdanger = 3.24 
made the formula for the solution change as t advanced though  
        tdanger = 3.24

T.3.b.i) When you're in the danger zone, your intuition may be thrown off

Here's another innocent-looking diffeq:

\(Clear[diffeq, y, t, f, starter]; 
f[t_, y_]\  = \ \ y\^\(2/3\); 


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

Go after a formula for the solution by separating variables:

Clear[yformula]; newleft=((y'[t])/(y[t]^(2/3)));newright = 1; neweqn = 
      ((∫_0^t) newleft \[DifferentialD]t /. y[0] -> 0) == 
        (∫_0^t) newright \[DifferentialD]t; \(ysol = Solve[neweqn, y[t]]; 


    formulay[t_] = 
      y[t] /. ysol[[1]]\)

Now plot:

endtime = 3; solutionplot = 
      Plot[formulay[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        AxesLabel -> {"t", "y[t]"}, 
        AspectRatio -> (1/GoldenRatio)];

And you check by showing this solution plot with a flow plot:

scalefactor = 0.4; {ylow, yhigh} = {0, 1}; \({tlow, thigh} = {0, endtime}; 


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


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

The solution is going with the flow.
Is there anything wrong here?

Only one thing:
There is another solution with the same starter data:
This is the silly solution 
          y[t] = 0 for all t ≥  0. 
This y[t] certainly solves the differential equation:

diffeq/.{y[t] -> 0,y'[t]->0,y[0] -> 0}

So the solution  y[t] = 0 deserves to be plotted too:

\(Clear[altsol]; 
altsol[t_] = 0; 


    altsolplot = 
      Plot[altsol[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], Red}}, 
        DisplayFunction -> Identity]; 


    bummer = Show[altsolplot, solutionplot, solutionflowplot, 
        Axes -> True, AxesLabel -> {"t", "y[t]"}, 
        AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> $DisplayFunction]; \)

You can see that the alternate silly solution is also going with the horizontal flow along the  x-axis.

T.b.ii) The danger zone strikes again

Take another look at the plot of the two solutions in part i) above:

Show[bummer];

When you throw your cork into this flow at {0, 0}, the cork cannot float away in two directions simultaneously.
Doesn't this blow the hell out of thinking of solution plots as paths of floating corks?

Yes it does. But before you give up on the idea of thinking of solution plots as paths of floating corks, check the danger zone for this differential equation:

\(Clear[diffeq, y, t, f, r, starter]; 
f[t_, y_]\  = y\^\(2/3\); 


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

gradf[t_, y_] = {(∂_t) f[t, y], (∂_y) f[t, y]}

There is a big fat singularity in gradf[x, y] at  y = 0,
So the danger zone is the line y[t] = 0:

ydangerous = 0; dangerzoneplot = 
      Graphics[{Thickness[0.005], Yellow, 
          Line[{{0, ydangerous}, {endtime, ydangerous}}]}]; \(dangerlabel = Graphics[Text["DANGER", {0.7\ endtime, ydangerous}]]; 
    

Show[bummer, dangerzoneplot, dangerlabel, 
      PlotLabel -> "Flow with Danger Zone"]; \)

The reason that your flow intuition doesn't quite work for this diffeq is that the cork was thrown into the flow in the danger zone.  In fact, the alternate solution stays in the danger zone all the time. The other one hustles its buns away from the danger zone.

T.3.b.iii)

Does this mean that when you are going with a differential equation with no danger zone, then you are free to use your intuition thinking of solution plots as the paths of floating corks?

Yes it does.
Fancy diffeq folks say that when there is no danger zone, solutions exist and are unique.