DE.06 Systems and Flows
 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)  Flows and their trajectories as pairs of solutions of a system of differential equations

B.1.a) Vector fields

A vector field is a function that spits out vectors. 
You make a 2D vector field by taking two regular functions, 
         m[x, y] and n[x, y], 
and throwing them into the two slots:

Clear[Field, m, n, x, y]; m[x_, y_] = 1.4 (y - 0.3); n[x_, y_] = 1.2 (x - 0.5); Field[x_, y_] = {m[x, y], n[x, y]}

At the sample point 
        {xsample,  ysample} = {0.3,  1.9}, 
the field vector is:

{xsample, ysample} = {0.3, 1.9}; Field[xsample, ysample]

You plot this vector with its tail at {xsample, ysample}:

ranger = 4; tailplot = 
      Graphics[{Red, PointSize[0.03], Point[{xsample, ysample}]}]; \(fieldvectorplot = 
      Arrow[Field[xsample, ysample], Tail -> {xsample, ysample}, 
        VectorColor -> Blue, HeadSize -> 0.8]; 


    Show[tailplot, fieldvectorplot, Axes -> True, 
      AxesLabel -> {"x", "y"}, 
      PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), ranger}}]; 
    \)

Here are some plots of some random field vectors for the given field
           Field[x, y] = {m[x, y], n[x, y]}:

ranger = 6; h = 3.0; {xsample, ysample} = {Random[Real, {-h, h}], 
        Random[Real, {-h, h}]}; tailplot = 
      Graphics[{Red, PointSize[0.03], Point[{xsample, ysample}]}]; \(fieldvectorplot = 
      Arrow[Field[xsample, ysample], Tail -> {xsample, ysample}, 
        VectorColor -> Blue, HeadSize -> 0.6]; 


    Show[tailplot, fieldvectorplot, Axes -> True, 
      AxesLabel -> {"x", "y"}, 
      PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), ranger}}]; 
    \)

B.1.a.i)

How do you see the whole thing?

You can plot a big sample of the vector field by plotting the vector, Field[x, y], with its tail at {x, y} for a selection of points {x, y}:

{xlow, xhigh} = {-3.0, 3.0}; {ylow, yhigh} = {-3.0, 3.0}; jump = 0.5; \(flowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> 0.5, HeadSize -> 0.4], {x, xlow, xhigh, 
          jump}, {y, ylow, yhigh, jump}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

A mad rush to the lower left and the upper right.

If you want to see unit vectors only, you can go with the option ScaleFactor->Normalize:

{xlow, xhigh} = {-3.0, 3.0}; {ylow, yhigh} = {-3.0, 3.0}; jump = 0.5; \(unitflowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> Normalize, HeadSize -> 0.4], {x, xlow, 
          xhigh, jump}, {y, ylow, yhigh, jump}]; 


    Show[unitflowplot, Axes -> True, 
      AxesLabel -> {"x", "y"}]; \)

You are free to take your pick.

B.1.a.ii) Floating corks

Stay with the flow plot:

Show[flowplot, Axes -> True, 
      AxesLabel -> {"x", "y"}];

What's a good way of interpreting the flow plot?

Look at what went into the plot.
The plot shows field vectors 
         Field[x, y] = {m[x, y], n[x, y]} 
plotted with their tails at {x, y} for a selection of points {x, y}.  The selection of points was:

{xlow, xhigh} = {-3.0, 3.0}; {ylow, yhigh} = {-3.0, 3.0}; jump = 0.5; \(pointplot = 
      Table[Graphics[{Red, PointSize[0.015], Point[{x, y}]}], {x, xlow, 
          xhigh, jump}, {y, ylow, yhigh, jump}]; 


    Show[pointplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

The plot of the whole flow shows field vectors, Field[x, y], plotted with their tails at the points {x, y} shown above.
Take a look:

Show[pointplot, flowplot, Axes -> True, 
      AxesLabel -> {"x", "y"}];

Now you're ready to interpret the plot.

Think of the whole xy-plane as fluid flowing with currents and eddies.
The vectors represent the flow of the fluid at their tails.
The length of each vector indicates the speed of the fluid flow at its tail.
The direction of each vector indicates the direction of the flow.
In the plot above, lots of the fluid is flowing off to the lower left and lots of the fluid is flowing off to the upper right.
The upshot:
The field vector Field[x, y] with tail at {x, y} represents the speed and the direction with which a cork at {x, y} moves away from {x, y} as it is caught by the flow.

B.1.a.iii)  Trajectories in the flow

Stay with the same vector field as above and look at this plot:

{xstarter, ystarter} = {2.0, -1.0}; endtime = 3.2; starterpoint = {xstarter, ystarter}; Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; \(ystartereqn = y[0] == ystarter; 


    approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    setup = Show[flowplot, starterplot, trajectoryplot, 
        PlotRange -> All, Axes -> True, 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction]; \)

Some folks like to call the path you see a "trajectory" of the vector field. Other folks call this path a "streamline" of the vector field. 
Just what is a trajectory in the field?

Take another look:

Show[setup];

This trajectory is part of the path of a cork thrown into the flow at the red dot. 
See some more:

{xstarter, ystarter} = {Random[Real, {1.5, 3}], 
        Random[Real, {-2, -1}]}; endtime = 3; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    \(trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; \)\)\(starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
      PlotRange -> {{xlow - 1, xhigh + 1}, {ylow - 1, yhigh + 1}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

And some more:

{xstarter, ystarter} = {Random[Real, {-2, -3}], 
        Random[Real, {2, 3}]}; endtime = 7; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    \(trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; \)\)\(starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
      PlotRange -> {{xlow - 1, xhigh + 1}, {ylow - 1, yhigh + 1}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Corks thrown into the flow at different points move on different paths.

B.1.a.iv)

Here is the same vector field:

Clear[Field, m, n, x, y]; m[x_, y_] = 1.4 (y - 0.3); \(n[x_, y_] = 1.2\ \((x - 0.5)\); 


    Field[x_, y_] = {m[x, y], n[x, y]}\)

How does this field govern the the shape of the trajectory of the cork dropped into the flow?

Answer:

If the cork finds itself at a position {x, y}, then it has to go with the flow.  The flow at {x, y} is in the direction of
        Field[x, y] = {m[x, y],  n[x, y]}:

Field[x, y]

{m[x, y], n[x, y]}

To go with the flow, the cork leaves {x, y} in the direction of the vector
        Field[x, y] = {m[x, y],  n[x, y]}.  
As the cork progresses, it corrects its direction instantaneously at each point to keep on the course that the flow determines.
As you saw above, the result is that if {x, y} is any point on the actual path of the cork, then the vector 
        Field[x, y] = {m[x, y],  n[x, y]}
with its tail at {x, y} is tangent to the actual path of the cork.
In this way, the vector field guides the cork as it moves along its path.

B.1.a.v) Pairs of solutions of a system of diffeq's as parameterizations of the trajectories

Here is the same vector field:

Clear[Field, m, n, x, y]; m[x_, y_] = 1.4 (y - 0.3); \(n[x_, y_] = 1.2\ \((x - 0.5)\); 


    Field[x_, y_] = {m[x, y], n[x, y]}\)

Use the given functions m[x, y] and n[x, y] to make this system of diffeq's:

\(\(
diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)\)

How are pairs of solutions {x[t],  y[t]} of this system of diffeq's related to the trajectories in the flow?

Answer:

Solutions {x[t], y[t]} of this system are parameterizations of the trajectories in the flow.
Reason: If at time t, the cork is at {x[t], y[t]}, then to go with the flow, the cork leaves {x[t], y[t]} in the direction of 
        Field[x[t], y[t]] = {m[x[t], y[t]],  n [x[t], y[t]}
at the speed indicated by the length of 
        {m[x[t], y[t]],  n [x[t], y[t]}.
        
On the other hand: 
If at time t, the cork is at {x[t],  y[t]}, then it leaves {x[t], y[t]} in the direction of 
        {x'[t],  y'[t]}
at the speed indicated by the length of 
         {x'[t],  y'[t]}. 
This tells you that
         {x'[t], y'[t]}  = 
         {m[x[t], y[t]],  n [x[t], y[t]}.
This is the same as:

\(\(
diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)\)

Different starter data on {x[0], y[0]} lead to different trajectories.

B.2)  Flow analysis of behavior of the solutions of the system 
                 x'[t]  =  m[x[t], y[t]]
                    y'[t]  =  n[x[t], y[t]]

B.2.a.i)

You enter a system of diffeq's by specifying two functions 
        m[x, y] and n[x, y]
like this:

\(Clear[m, n, x, y]; 
m[x_, y_] = \(-0.82\)\ x\  - 0.17\ y\)n[x_, y_] = 0.50 x  - 0.29 y

And then you use these functions to make your diffeq system
        x'[t]=m[x[t], y[t]]
        y'[t]=n[x[t], y[t]]
with specified starting values for x[0] and y[0];

\(diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)

Use flow plots to estimate what happens to the solutions {x[t], y[t]} as t grows and grows.

Go with a scaled plot of Field[x, y] = {m[x, y], n[x, y]}:

\(Clear[m, n, x, y]; 
m[x_, y_] = \(-0.82\)\ x - 0.17\ y; 

    n[x_, y_] = 0.50\ x - 0.29\ y; 

Clear[Field]; 

    Field[x_, y_] = {m[x, y], n[x, y]}; 
{xlow, xhigh} = {\(-4\), 4}; 

    \({ylow, yhigh} = {\(-4\), 4}; \)\)\(flowplot = 
      Table[Arrow[{m[x, y], n[x, y]}, Tail -> {x, y}, 
          VectorColor -> Blue, ScaleFactor -> 0.35, 
          HeadSize -> 0.25], {x, xlow, xhigh, 0.50}, {y, ylow, yhigh, 
          0.50}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

All the trajectories seem to be swirling in toward {0, 0}.
Check this out by dropping in a cork at a random point;

{xstarter, ystarter} = {Random[Real, {xlow, xhigh}], 
        Random[Real, {ylow, yhigh}]}; endtime = 10; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    \(trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; \)\)\(starterplot = Graphics[{Red, PointSize[0.05], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
      PlotRange -> {{xlow - 1, xhigh + 1}, {ylow - 1, yhigh + 1}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Yep, the trajectories wind their way to {0, 0}.
See three more trajectories:

{xstarter1, ystarter1} = {Random[Real, {xlow, xhigh}], 
        Random[Real, {ylow, 0}]}; {xstarter2, ystarter2} = {Random[Real, {xlow, 0}], 
        Random[Real, {0, yhigh}]}; {xstarter3, ystarter3} = {Random[Real, {0, xhigh}], 
        Random[Real, {ylow, 0}]}; endtime = 8; starterpoint1 = {xstarter1, ystarter1}; starterpoint2 = {xstarter2, ystarter2}; starterpoint3 = {xstarter3, ystarter3};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn1 = x[0] == xstarter1; ystartereqn1 = y[0] == ystarter1; xstartereqn2 = x[0] == xstarter2; ystartereqn2 = y[0] == ystarter2; xstartereqn3 = x[0] == xstarter3; ystartereqn3 = y[0] == ystarter3; ndssol1 = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn1, ystartereqn1}, {x[t], y[t]}, {
          t, 0, endtime}]; ndssol2 = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn2, ystartereqn2}, {x[t], y[t]}, {
          t, 0, endtime}]; \(ndssol3 = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn3, ystartereqn3}, {x[t], y[t]}, {
          t, 0, endtime}]; 

Clear[x1, x2, x3, y1, y2, y3]; 
    
{x1[t_], y1[t_]} = {
        x[t] /. ndssol1[[1]], 
        y[t] /. ndssol1[[1]]}; 

    \({x2[t_], y2[t_]} = {
        x[t] /. ndssol2[[1]], 
        y[t] /. ndssol2[[1]]}; \)\)\({x3[t_], y3[t_]} = {
        x[t] /. ndssol3[[1]], 
        y[t] /. ndssol3[[1]]}; 


    Clear[trajectory1, trajectory2]; 
trajectory1[t_] = {x1[t], y1[t]}; 

    \(trajectory2[t_] = {x2[t], y2[t]}; \)\)trajectory3[t_] = {x3[t], y3[t]}; \(trajectoryplots = 
      ParametricPlot[{trajectory1[t], trajectory2[t], trajectory3[t]}, {t, 0, 
          endtime}, PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplots = {Graphics[{Red, PointSize[0.05], Point[starterpoint1]}], 
        Graphics[{Red, PointSize[0.05], Point[starterpoint2]}], 
        Graphics[{Red, PointSize[0.05], Point[starterpoint3]}]}; 


    Show[flowplot, starterplots, trajectoryplots, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

The trajectories are sucked toward {0, 0}. 
See corks flow on the last trajectories:

\(Clear[trajplots, s]; 

    trajplots[s_] := 
      ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}}, {t, 0, 
          s}, PlotStyle -> {{Thickness[0.01], CadmiumOrange}}, 
        DisplayFunction -> Identity]; 


    Table[Show[flowplot, trajplots[s], starterplots, Axes -> True, 
        PlotLabel -> s\ "= endtime", 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction], {s, 0.5, endtime, 
        \(endtime - 0.5\)\/4}]; \)

The corks are sucked into {0, 0}.   
The visual evidence is that when you go with:

\(diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)

then when t gets big, both x[t] and y[t] get close to 0.

B.2.a.ii)

Stay with the same setup as in part i) immediately above and look at this trajectory shown with the flow plot:

{xstarter, ystarter} = {3.5, -1}; endtime = 10; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(ndssol = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 


    Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.05], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
      PlotRange -> {{xlow - 1, xhigh + 1}, {ylow - 1, yhigh + 1}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

This is a plot of the trajectory starting at the indicated point. 
Plot the individual solutions x[t] and y[t] corresponding to the given starter data.

The solutions come in pairs {x[t], y[t]} for this trajectory. 
Here are the individual solution plots:

xplot = 
      Plot[x[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Magenta}}, 
        AxesLabel -> {"t", "x[t]"}, PlotRange -> All, 
        AspectRatio -> (1/GoldenRatio), 
        PlotLabel -> "x solution", 
        Epilog -> {Red, PointSize[0.05], Point[{0, xstarter}]}];

yplot  = 
      Plot[y[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Blue}}, 
        AxesLabel -> {"t", "y[t]"}, PlotRange -> All, 
        AspectRatio -> (1/GoldenRatio), 
        PlotLabel -> "y solution", 
        Epilog -> {Red, PointSize[0.05], Point[{0, ystarter}]}];

See them side-by-side:

Show[GraphicsArray[{xplot, yplot}]];

The corresponding trajectory {x[t], y[t]} plots out this way:

trajectoryplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], CadmiumOrange}}, 
        AxesLabel -> {"x", "y"}, PlotRange -> All, 
        PlotLabel -> "Trajectory plot", 
        Epilog -> {Red, PointSize[0.05], Point[{x[0], y[0]}]}];

Here are plots of x[t], and y[t], and of the trajectory, {x[t], y[t]} all at once:

Show[GraphicsArray[{{xplot, yplot}, {trajectoryplot}}]];

The pair {x[t], y[t]} IS the solution.
The trajectory IS NOT the solution, but it does do a good job of showing off what the solutions do.
You can see from the trajectory plot that the corresponding x[t] solution plot goes down as time t advances. 
You can also see from the trajectory plot that the corresponding y[t] solution plot goes up and then down as time t advances.

B.2.a.iii)

Keep everything the same as immediately above.
Can you see the solution plots and the trajectory all in one plot?

If you are willing to look at a 3D plot, then you can.
Here you go:

threeDstarterpoint = {0, 3.5, -1}; threeDstarterplot = 
      Graphics3D[{Red, PointSize[0.03], Point[threeDstarterpoint]}]; allthree = 
      ParametricPlot3D[{t, x[t], y[t]}, {t, 0, endtime}, 
        AxesLabel -> {"t", "x[t]", "y[t]"}, 
        DisplayFunction -> Identity]; \(taxis = Graphics3D[{Green, Line[{{0, 0, 0}, {endtime, 0, 0}}]}]; 

    xaxis = Graphics3D[{Green, Line[{{0, 0, 0}, {0, 3.5, 0}}]}]; 

    yaxis = Graphics3D[{Green, Line[{{0, 0, \(-1\)}, {0, 0, 1}}]}]; 


    everything = 
      Show[allthree, threeDstarterplot, taxis, xaxis, yaxis, 
        Boxed -> False, PlotRange -> All, ViewPoint -> CMView, 
        DisplayFunction -> $DisplayFunction]; \)

To see the trajectory plot, look at this 3D curve from a view point down the t - axis:

Show[everything, ViewPoint -> {25, 0, 0}];

Compare:

Show[trajectoryplot];

To see the x[t] solution plot, look at the 3D curve from a view point down the y - axis:

Show[everything, ViewPoint -> {0, 0, 25}];

Compare:

Show[xplot];

To see the y[t] solution plot, look at the 3  D curve from a view point down the negative x- axis:

Show[everything, ViewPoint -> {0, -25, 0}];

Compare:

Show[yplot];

B.2.b.i) Oscillations

Here's a new system of differential equations
         x'[t]=1.8 y[t]
         y'[t]=-3.2 x[t]
        with given starter data on x[0] and y[0].
Do a flow analysis to predict what happens to the solutions {x[t], y[t]} as t grows and grows.

Look at the system:
          x'[t]=1.8 y[t]
          y'[t]=-3.2 x[t].
Enter the corresponding functions m[x, y] and n[x, y] this way:

Clear[m, n, x, y, t, starterx, startery];
m[x_, y_] = 1.8 y;
n[x_, y_] = -3.2 x;

diffeqsystem={x'[t],y'[t]}=={m[x[t], y[t]], n[x[t], y[t]]};
({{diffeqsystem}})

Now that you have entered m[x, y] and n[x, y], you can look at a scaled plot of Field[x, y] = {m[x, y], n[x, y]}:

\(Clear[Field]; 
Field[x_, y_] = {m[x, y], n[x, y]}\){xlow, xhigh} = {-4, 4}; {ylow, yhigh} = {-5, 5}; jump = 1; \(flowplot = 
      Table[Arrow[{m[x, y], n[x, y]}, Tail -> {x, y}, 
          VectorColor -> Blue, ScaleFactor -> 0.25, 
          HeadSize -> 0.5], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

All the trajectories seem to be swirling around {0, 0}.
Check this out by dropping in a cork at {3, -2};

endtime = 4; {xstarter, ystarter} = {3, -2}; endtime = 10; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(ndssol = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 


    Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.05], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
      PlotRange -> {{xlow - 1, xhigh + 1}, {ylow - 1, yhigh + 1}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Yep, the trajectory swirls around {0, 0}
See three more trajectories starting at some partially random points:

{xstarter1, ystarter1} = {Random[Real, {-6, -2}], 0.0}; {xstarter2, ystarter2} = {-1.7, Random[Real, {0, 6}]}; {xstarter3, ystarter3} = {0.0, Random[Real, {1, 6}]}; endtime = 5; starterpoint1 = {xstarter1, ystarter1}; starterpoint2 = {xstarter2, ystarter2}; starterpoint3 = {xstarter3, ystarter3};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn1 = x[0] == xstarter1; ystartereqn1 = y[0] == ystarter1; xstartereqn2 = x[0] == xstarter2; ystartereqn2 = y[0] == ystarter2; xstartereqn3 = x[0] == xstarter3; ystartereqn3 = y[0] == ystarter3; approxsolutions1 = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn1, ystartereqn1}, {x[t], y[t]}, {
          t, 0, endtime}]; approxsolutions2 = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn2, ystartereqn2}, {x[t], y[t]}, {
          t, 0, endtime}]; \(approxsolutions3 = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn3, ystartereqn3}, {x[t], y[t]}, {
          t, 0, endtime}]; 

Clear[x1, x2, x3, y1, y2, y3]; 
    
{x1[t_], y1[t_]} = {
        x[t] /. approxsolutions1[[1]], 
        y[t] /. approxsolutions1[[1]]}; 
    
\({x2[t_], y2[t_]} = {
        x[t] /. approxsolutions2[[1]], 
        y[t] /. approxsolutions2[[1]]}; 
    \)\)\({x3[t_], y3[t_]} = {
        x[t] /. approxsolutions3[[1]], 
        y[t] /. approxsolutions3[[1]]}; 
    

Clear[trajectory1, trajectory2]; 
trajectory1[t_] = {x1[t], y1[t]}; 
    
\(trajectory2[t_] = {x2[t], y2[t]}; \)\)trajectory3[t_] = {x3[t], y3[t]}; trajectoryplots = 
      ParametricPlot[{trajectory1[t], trajectory2[t], trajectory3[t]}, {t, 0, 
          endtime}, PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; \(starterplots = {Graphics[{Red, PointSize[0.06], Point[starterpoint1]}], 
        Graphics[{Red, PointSize[0.06], Point[starterpoint2]}], 
        Graphics[{Red, PointSize[0.06], Point[starterpoint3]}]}; 


    Show[flowplot, starterplots, trajectoryplots, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

The trajectories swirl around {0, 0}. 
See them move as time advances:

\(Clear[trajplots, s]; 

    trajplots[s_] := 
      ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}}, {t, 0, 
          s}, PlotStyle -> {{Thickness[0.015], CadmiumOrange}}, 
        DisplayFunction -> Identity]; 


    Table[Show[flowplot, starterplots, trajplots[s], Axes -> True, 
        PlotLabel -> s\ "= time", 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction], {s, 0.5, endtime, 
        \(endtime - 0.5\)\/4}]; \)

The corks swirl around {0, 0}. 
Once it finishes one cycle, it repeats its trip over and over. 
Play with some other starting points. 

The visual evidence is that when you go with:

Clear[m, n, x, y, t, starterx, startery];
m[x_, y_] = 1.8 y;
n[x_, y_] = -3.2 x;

diffeqsystem={x'[t],y'[t]}=={m[x[t], y[t]], n[x[t], y[t]]};
({{diffeqsystem}})

Then when t gets bigger and bigger, all trajectories swirl around {0, 0} over and over. 
This means that plots of the individual solutions x[t] and y[t] oscillate like sine waves.
Here's the plot of the solution x[t] corresponding to some random starting data on {x[0], y[0]} followed by the plot of the companion solution y[t]:

endtime = 8; {xstarter, ystarter} = {Random[Real, {xlow, xhigh}], 
        Random[Real, {ylow, yhigh}]}; endtime = 10; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
xplot = 
      Plot[x[t], {t, 0, endtime}, 
        PlotStyle -> {{Magenta, Thickness[0.01]}}, 
        AxesLabel -> {"t", "x[t]"}, 
        AspectRatio -> 1\/GoldenRatio, 
        PlotLabel -> "x[t] solution plot"]; 


    yplot = Plot[y[t], {t, 0, endtime}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"t", "y[t]"}, 
        AspectRatio -> 1\/GoldenRatio, 
        PlotLabel -> "y[t] solution plot"]; \)

Here's the trajectory {x[t], y[t]} with its starter point:

\(trajplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], CadmiumOrange}}, 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> Identity]; 


    trajectory = 
      Show[trajplot, Graphics[{PointSize[0.08], Red, Point[starterpoint]}], 
        DisplayFunction -> $DisplayFunction]; \)

See them all together:

Show[GraphicsArray[{{xplot, yplot}, {trajectory}}]];

B.2.b.ii)

Stay with the same diffeq system as in part i):

\(Clear[m, n, x, y, t, starterx, startery]; 
m[x_, y_]\  = \ 1.8\ y; 

    n[x_, y_]\  = \ \(-3.2\)\ x; 


    diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)

If you want a plot of only y[t], then can you get by with given starter data on y[0] without any given starter data on x[0]?

You cannot get by with given starter data on y[0] without any given starter data on x[0].
Reason:
The solution plots come from the trajectories and in order to get a trajectory you have to know where the cork was thrown into the flow.
To have any chance of doing this, you need starter data on both x[0] and y[0].

B.3)  Flow analysis of  the unforced linearoscillator differential equation by converting it to a system of two first order differential equations

B.3.a.i) The undamped linear oscillator differential equation

Here's a prototype of the undamped linear oscillator differential equation:

\(Clear[t, y]; 

    oscillatordiffeq\  = \ \((\(\(y'\)'\)[t]\  + \ 1.2\ y[t]\  == \ 0)\)\)

See a sample solution:

Clear[y1];
endtime = 17;
y1[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == 4,y'[0]==-2},y[t],{t, 0, endtime}][[1]];

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

It oscillates just as a sine wave oscillates.
See three more sample solutions:

Clear[y1, y2, y3, y, t];
y1[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == -3,y'[0]==-2},y[t],{t, 0, endtime}][[1]];
y2[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == 0,y'[0]==1},y[t],{t, 0, endtime}][[1]];
y3[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == -2,y'[0]==2},y[t],{t, 0, endtime}][[1]];

solutionplots = 
        Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
          PlotStyle 
            -> {{Thickness[0.01], Blue}, {Thickness[0.01], Red}, {
                Thickness[0.01], Orange}}, 
          AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
          AxesLabel -> {"t", ""}];

They all oscillate just the way a sine wave oscillates.

B.3.a.ii) Starter data on y[0] and y'[0]

Go back and look at the code that produced these solutions and notice that starter data on y[0] and y'[0] are given for each solution. Why is this natural?

To see why it is as natural as apple pie, you just make this oscillator differential equation into a system of two first order differential equations. 
To do this you make up a new function.  Call it x[t], and set 
             x[t]=y'[t]. 
And then:
→ Every time you see a y'[t] in the differential equation, replace it with an x[t].  
→ Every time you see a y''[t] in the equation, replace it with an x'[t].  
Take a look.

\(Clear[x]; 

    ColumnForm[
      Thread[{new\  = \ \((\(y'\)[t]\  == \ x[t])\), 

          oscillatordiffeq /. {\(y'\)[t] -> x[t], 
              \(\(y'\)'\)[t] -> \(x'\)[t]}}]\ ]\)

Clean this up by putting the x'[t] on the left and the y[t] on the right:

\(\(
Clear[m, n, x, y, t, starterx, startery]; 

    m[x_, y_]\  = \ \(-\ 1.2\)\ y; 
n[x_, y_]\  = \ \ x; 


    \((diffeqsystem\  = 
          \ {\(x'\)[t]\  == \ m[x[t], y[t]], \(y'\)[t]\  == \ n[x[t], y[t]]})
        \) // ColumnForm\)\)

This system is exactly the same as:

oscillatordiffeq

Any y[t] that solves the diffeq system also solves the oscillator diffeq. 
To nail down a y[t] for the system, you need starter data on x[0] and y[0]. 
But because 
          x[t]=y'[t], 
the starter data on x[0] is actually starter data on y'[0]. 
That's why you need starter data on both y[0] and y'[0] to nail down a solution of the oscillator diffeq.

Play by going with random starter data on y[0] and y'[0] and re-running several times:

\({startery, starteryprime}\  = 
      \ {Random[Real, {\(-2\), 2}], Random[Real, {\(-2\), 2}]}; 
Clear[y1]; 
    
endtime\  = \ 17; 


    y1[t_]\  = \ 
      y[t] /. \(NDSolve[
            
{oscillatordiffeq, y[0]\  == \ startery, \ 
              \(y'\)[0]\  == \ starteryprime}, 
\ \ \ \ \ \ \ \ \ 
            y[t], {t, 0, endtime}]\)[\([1]\)]; 


    Plot[y1[t], {t, 0, endtime}, PlotStyle -> 
	{{Thickness[0.01], Red}}, 
      PlotRange -> {\(-3\), 3}, 
	AspectRatio -> 1\/GoldenRatio, 
      PlotRange -> All, AxesLabel -> {"t", ""}]; \)

B.3.a.iii) Converting the unforced undamped oscillator diffeq to a system to see why solutions must oscillate like sine waves

Here's another undamped, unforced oscillator diffeq:

\(Clear[t, y]; 

    oscillatordiffeq\  = \ \((\(\(y'\)'\)[t]\  + \ 2.2\ y[t]\  == \ 0)\)\)

Explain how you can use a flow plot to tell that all solutions of this diffeq oscillate like sine waves.

Convert it to a system and look at a flow plot:

\(Clear[x]; 

    ColumnForm[
      Thread[{new\  = \ \((\(y'\)[t]\  == \ x[t])\), 

          oscillatordiffeq /. {\(y'\)[t] -> x[t], 
              \(\(y'\)'\)[t] -> \(x'\)[t]}}]\ ]\)

Clean this up:

\(Clear[m, n, x, y, t, starterx, startery]; 

    m[x_, y_]\  = \ \(-\ 2.2\)\ y; 
n[x_, y_]\  = \ \ x; 


    \((diffeqsystem\  = 
          \ {\(x'\)[t]\  == \ m[x[t], y[t]], \(y'\)[t]\  == \ n[x[t], y[t]]})
        \) // ColumnForm\)

See the flow:

\(Clear[Field]; 
Field[x_, y_] = {m[x, y], n[x, y]}\){xlow, xhigh} = {-4, 4}; {ylow, yhigh} = {-4, 4}; jump = 0.8; \(flowplot = 
      Table[Arrow[{m[x, y], n[x, y]}, Tail -> {x, y}, 
          VectorColor -> Blue, ScaleFactor -> 0.4, 
          HeadSize -> 0.5], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

Your eyes tell you that the trajectories swirl around {0, 0}. On a specific trajectory, y[t] goes up and down over and over. That's why the oscillator oscillates. 
See a random trajectory:

endtime = 4; {xstarter, ystarter} = {Random[Real, {xlow, xhigh}], 
        Random[Real, {ylow, yhigh}]}; endtime = 10; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(ndssol = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 


    Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
      PlotRange 
        -> {{xlow - 1.5, xhigh + 1.5}, {ylow - 1.5, yhigh + 1.5}}, 
      Axes -> True, AxesLabel -> {"x = y'", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

See a solution y[t] of the linear oscillator diffeq oscillate up and down. 
Check it out:

yfromsystem = 
      Plot[y[t], {t, 0, endtime}, PlotStyle -> {{Red, Thickness[0.02]}}, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "From the system"];

Compare with a direct NDSolve:

Clear[y1, y];
y1[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == ystarter,y'[0]==xstarter},y[t],{t, 0, endtime}][[1]];

\(directplot = 
        Plot[y1[t], {t, 0, endtime}, 
		
          PlotStyle -> {{Thickness[0.01], Blue}}, 
          AspectRatio -> 1\/GoldenRatio, PlotRange -> All, 
          AxesLabel -> {"t", ""}, 
          PlotLabel -> "From NDSolve"]\);

See them both:

Show[yfromsystem, directplot, PlotLabel -> "Both"];

Dead ringers. And this is no accident. It will happen everytime.

Check out other solutions corresponding to random starting data on y[0] and y'[0]:

{startery, starteryprime} = {Random[Real, {xlow, xhigh}], 
          Random[Real, {ylow, yhigh}]};

Clear[y, t, approxy];
endtime = Random[Real, {10, 25}];

Clear[y1];
y1[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == startery,y'[0]==starteryprime},y[t],{t, 0, endtime}][[1]];

\(Plot[y1[t], {t, 0, endtime}, 
	
        PlotStyle -> {{Thickness[0.015], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}]\);

Rerun until you feel good.

B.3.b) The damped unforced oscillator

You get a damped unforced oscillator by inserting a b y'[t] term (with b > 0) into the oscillator diffeq:

\(Clear[t, y]; 
b\  = \ 1.8; 

    oscillatordiffeq\  = \ 
      \((\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ 6.4\ y[t]\  == \ 0)\)\)

Take a look at some sample solutions coming from different starter data on y[0] and y'[0]:

Clear[y1, y2, y3, y, t];
endtime = 6;
y1[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == 3,y'[0]==-2},y[t],{t, 0, endtime}][[1]];y2[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == 0,y'[0]==1},y[t],{t, 0, endtime}][[1]];y3[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == -2,y'[0]==2},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 -> All, 
          AxesLabel -> {"t", ""}];

These sample solutions look like damped (squished) sine waves. They try to oscillate but are squashed out as time goes on.
Convert this damped oscillator diffeq into a diffeq system and use the resulting flow plot of the system to explain the behavior you see above.

Convert it to a system by replacing all y'[t]'s with x[t] and all (y^′)′[t]'s  by x'[t]:

\(Clear[x]; 

    ColumnForm[
      Thread[{new\  = \ \((\(y'\)[t]\  == \ x[t])\), 

          oscillatordiffeq /. {\(y'\)[t] -> x[t], 
              \(\(y'\)'\)[t] -> \(x'\)[t]}}]\ ]\)

Enter the correct m[x, y] and n[x, y] and clean this up:

\(Clear[m, n, x, y, t, starterx, startery]; 

    m[x_, y_]\  = \ \(-\ 6.4\)\ y\  - \ 1.9\ x; 
n[x_, y_]\  = \ \ x; 


    \((diffeqsystem\  = 
          \ {\(x'\)[t]\  == \ m[x[t], y[t]], \(y'\)[t]\  == \ n[x[t], y[t]]})
        \) // ColumnForm\)

Look at the flow plot:

\(Clear[Field]; 
Field[x_, y_] = {m[x, y], n[x, y]}\){xlow, xhigh} = {-3, 3}; {ylow, yhigh} = {-3, 3}; jump = ((xhigh - xlow)/12); \(flowplot = 
      Table[Arrow[{m[x, y], n[x, y]}, Tail -> {x, y}, 
          VectorColor -> Blue, ScaleFactor -> 0.15, 
          HeadSize -> 0.4], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

Your eyes tell you that the trajectories swirl TOWARD {0, 0}. 

On a specific trajectory, y[t] goes up and down over and over but it gets closer and closer to 0 as time goes on. 

That's why the damped oscillator oscillates while it is being squashed. See a random trajectory:

endtime = 9; {xstarter, ystarter} = {Random[Real, {xlow, xhigh}], 
        Random[Real, {ylow, yhigh}]}; endtime = 10; starterpoint = {xstarter, ystarter};

Clear[x, y, t];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(ndssol = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 


    Clear[trajectory]; 
trajectory[t_] = {x[t], y[t]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[flowplot, starterplot, trajectoryplot, 
	
      PlotRange -> {{xlow - 1, xhigh + 1}, {ylow - 1, yhigh + 1}}, 
      Axes -> True, AxesLabel -> {"x = y'", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

See the y[t] coming from the last plot oscillate up and down on its way to 0:

yplot = 
      Plot[y[t], {t, 0, endtime}, PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"t", "y[t]"}, PlotRange -> All, 
        AspectRatio -> (1/GoldenRatio), 
        PlotLabel -> "From the system"];

Check out other solutions of

\(Clear[y]; 
	oscillatordiffeq\)

corresponding to random starting data on y[0] and y'[0]:

{startery, starteryprime} = {Random[Real, {-3, 3}], 
          Random[Real, {-3, 3}]};

Clear[y, t, approxy];
endtime = 8;
Clear[y1, y];
y1[t_]=y[t]/.NDSolve[{oscillatordiffeq,y[0] == startery,y'[0]==starteryprime},y[t],{t, 0, endtime}][[1]];

Plot[y1[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}];

Rerun until you are in the swing.

B.4) The danger zone: Where your cork flow intuition can be blown to bits

B.4.a) Weirdness in the Danger Zone

Here's a little system:

\(Clear[m, n, x, y, t]; 
m[x_, y_]\  = \ x\^\(2/3\); 

    n[x_, y_]\  = \ y\^\(2/3\)\ ; 


    diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)

And a flow plot:

\(Clear[Field]; 
Field[x_, y_] = {m[x, y], n[x, y]}; 

    \(scalefactor = 0.2; \)\){xlow, xhigh} = {0, 5}; {ylow, yhigh} = {0, 5}; \(flowplot = 
      Table[Arrow[{m[x, y], n[x, y]}, Tail -> {x, y}, 
          VectorColor -> Blue, ScaleFactor -> Normalize, 
          HeadSize -> 0.3], {x, xlow, xhigh, 0.5}, {y, ylow, yhigh, 
          0.5}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

Here are three formulas for pairs of solutions popping out of {0, 0}:

\(Clear[x1, y1, x2, y2, x3, y3, t]; 
    
{x1[t_], y1[t_]} = {t\^3\/27, t\^3\/27}\){x2[t_], y2[t_]} = {((t^3)/27), 0}{x3[t_], y3[t_]} = {0, ((t^3)/27)}

Short hand calculations show each pair solves:

diffeqsystem

Watch their trajectories go with the flow:

\(Clear[trajplots, endtime]; 

    trajplots[endtime_] := 
      ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}}, {t, 0, 
          endtime}, PlotStyle -> {{Thickness[0.015], CadmiumOrange}}, 
        DisplayFunction -> Identity]; 


    Table[Show[flowplot, trajplots[endtime], Axes -> True, 
        PlotLabel -> endtime\ "= endtime", 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction], {endtime, 2, 5, 1}]; \)

Wait a minute!
You drop the cork into the flow at {0, 0} and it floats off in three different directions simultaneously! 
Doesn't this shoot to hell all the intuition you get by thinking of trajectories as corks going with the flow?

Not if you are aware of the strange things that can and sometimes do happen in the danger zone.

To locate the danger zone, look again at the setup:

\(Clear[m, n, x, y, t]; 
m[x_, y_] = x\^\(2/3\); 

    n[x_, y_] = y\^\(2/3\); 


    DEsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[DEsystem]]\)

Take a look at the gradients of m[x, y] and n[x, y]:

\(Clear[gradm]; 

    gradm[x_, y_] = {∂\_x m[x, y], ∂\_y m[x, y]}\)

\(Clear[gradn]; 

    gradn[x_, y_] = {∂\_x n[x, y], ∂\_y n[x, y]}\)

The danger zone consists of the singularities (blow-ups) of gradm[x, y] or the singularities (blow-ups) of gradn[x, y].

By inspection, you can see that the danger zone consists of the lines
                     x = 0 and y = 0. 
 See the three trajectories plotted with the danger zone:

dangerzone = {
        Graphics[{Yellow, Thickness[0.02], Line[{{0, 0}, {5, 0}}]}], 
        Graphics[{Yellow, Thickness[0.02], Line[{{0, 0}, {0, 5}}]}]}; \(dangerlabel = {Graphics[Text["DANGER", {5\/2, 0}]], 
        Graphics[Text["  DANGER", {0, 5\/2}]]}; 


    Table[Show[flowplot, flowplot, dangerzone, dangerlabel, 
        trajplots[endtime], Axes -> True, 
        PlotLabel -> endtime\ "= endtime", 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction], {endtime, 2, 5, 1}]; \)

All three solutions start in the danger zone and two of them never leave the danger zone. 

Morals:
Whenever you start a trajectory in a danger zone you should be on the lookout for something that violates your intuition.
Anytime your trajectory enters a danger zone, expect the unexpected.

B.4.b)

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

Go with a cleared diffeq system:

\(Clear[diffeqsystem, x, y, m, n, t]; 

    DEsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[DEsystem]]\)

Because 
       x'[t]=m[x[t], y[t]]
and 
       y'[t]=n[x[t], y[t]]
the chain rule from vector calc guarantees that
       x''[t]=gradm[x[t], y[t]].{x'[t],y'[t]}.  
       y''[t]=gradn[x[t], y[t]].{x'[t],y'[t]}.

Now here's the key point: 
If a trajectory {x[t], y[t]} hits the danger zone (i.e. a singularity of gradm[x, y] or of gradn[x, y]), at a time t  =  tdanger, then the facts that 
       x''[tdanger]=gradm[x[tdanger], y[tdanger]].{x'[tdanger],y'[tdanger]}.  
       y''[tdanger]=gradn[x[tdanger], y[tdanger]].{x'[tdanger],y'[tdanger]}.
 opens up the possibility of singularities of 
        (x^′)′[t] or (y^′)′[t] at t = tdanger.
When this happens things can get out of hand! 
And they did in part i) above.