DE.04 Modern DiffEq Issues
 BASICS

Needs["Graphics`Colors`"];
Graphics`Colors`GosiaGreen=RGBColor[0, 0.392187, 0];
Needs["Graphics`FilledPlot`"];

If[MemberQ[{"Global`","Graphics`Common`GraphicsCommon`","Vector3D`"}, 
    Context[ScaleFactor]],Remove["*`ScaleFactor"]];
<<"Graphics`Legend`";
<<"Graphics`Common`GraphicsCommon`";
Remove["*`ScaleFactor"];

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) Euler's method of faking the plot of the solution of a differential equation
        and how it highlights the fundamental issue of diffeq

Here's an innocent-looking differential equation:
       y'[x]  =  y[x] (4 Cos[x] -  y[x])
       with y[0]  =  starter  =  1. 
If you want to get a plot of the solution, you ask Mathematica for a formula for the solution:

Clear[x, y];
starter = 1;
DSolve[{y'[x]==y[x] (4 Cos[x]^2 - y[x]),y[0] == starter},y[x],x]

This output is useless.  
Try as it might, Mathematica could not come up with a useful formula for the solution y[x].  But there is a way to fake the plot of the solution even when you can't get a formula for it.
It comes from this code which you can play with.
Don't  worry  about  the  individual  statements  in  the  code  right  now.

a = 0; b = 5; starter = 1.7; y[a] = starter; beginner = {a, y[a]}; \(Clear[x, y, next, jump, iterations, point, k]; 

    jump[iterations_] := \(b - a\)\/iterations; 

    \(next[{x_, y_}, iterations_] := {x, y} + 
        jump[iterations]\ {1, y\ \((4\ Cos[x]\^2 - y)\)}; \)\)point[0, iterations_] = beginner; point[k_, iterations_] := 
      point[k, iterations] = 
        N[next[point[k - 1, iterations], iterations]]; \(Clear[Euler]; 

    Euler[iterations_] := 
      Show[Graphics[{Thickness[0.01], Red, 
            Line[Table[point[k, iterations], {k, 0, iterations}]]}], 
        PlotRange -> All, Axes -> True, 
        AxesOrigin -> beginner, 
        AxesLabel -> {"x", "Euler fake"}, 
        PlotLabel -> iterations\ " iterations", 
        DisplayFunction -> Identity]; 

    \(primitive = Show[Euler[10], DisplayFunction -> $DisplayFunction]; 
    \)\)

A primitive stick figure resulting from 10 iterations. 
Here's what you get with 20 iterations:

better = Show[Euler[20], DisplayFunction -> $DisplayFunction];

50 iterations:

evenbetter = 
      Show[Euler[50], DisplayFunction -> $DisplayFunction];

100 iterations:

darngood = 
      Show[Euler[100], DisplayFunction -> $DisplayFunction];

Notice that the last two plots are almost the same.  
This is your signal that they are both reasonably good fake plots of the solution of the differential equation.

Show[GraphicsArray[{{primitive, better}, {evenbetter, darngood}}]];

The more iterations you use, the better the fake plot of the solution of the given differential equation.

B.1.a)

Here is the exponential differential equation 
         y'[x]  =  0.3  y[x] 
         with  y[0]  =  4 
and a formula for its solution:

thisr = 0.3; thisstarter = 4; Clear[x, y, r, starter];
solution=DSolve[{y'[x]==r y[x],y[0] == starter},y[x],x]/.{r -> thisr, starter -> thisstarter};
y[x_] = 
        y[x] /. solution[[1]]

and a plot for 0 ≤ x ≤ 8:

realplot = 
      Plot[y[x], {x, 0, 8}, PlotStyle -> {{Blue, Thickness[0.015]}}, 
        PlotRange -> All, AspectRatio -> (1/GoldenRatio), 
        AxesLabel -> {"x", ""}];

Compare the Euler fake plots of the solution with the true plot. Discuss what you see.

Copy, paste, and edit the code used above to get a Euler faker of a plot of the solution of the differential equation 
         y'[x]  =  0.3 y[x]
         with  y[0]  =  4 
and show the fake plots together with the plot of the true solution.

a = 0; b = 8; starter = 4; y[a] = starter; beginner = {a, y[a]}; \(Clear[x, y, next, jump, iterations, point, k]; 

    jump[iterations_] := \(b - a\)\/iterations; 

    next[{x_, y_}, iterations_] := {x, y} + jump[iterations]\ {1, 0.3\ y}; 

    
point[0, iterations_] = beginner; 

    point[k_, iterations_] := 
      \(point[k, iterations] = 
        N[next[point[k - 1, iterations], iterations]]\); 

Clear[Euler]; 

    Euler[iterations_] := 
      Show[Graphics[{Thickness[0.01], Red, 
            Line[Table[point[k, iterations], {k, 0, iterations}]]}], 
        PlotRange -> All, Axes -> True, 
        AxesOrigin -> beginner, AxesLabel -> {"x", ""}, 
        PlotLabel -> iterations\ " iterations", 
        DisplayFunction -> Identity]; 


    rough = Show[realplot, Euler[10], PlotLabel -> "10 iterations", 
        DisplayFunction -> $DisplayFunction]; \)

Step up the number of iterations to 20:

\(\(better\  = \ 
      Show[realplot, Euler[20], 
\ \ \ \ \ \ \ \ \ \ \ \ 
        PlotLabel -> "20 iterations", 
\ \ \ \ \ \ \ \ \ \ \ \ 
        DisplayFunction -> $DisplayFunction]; \)\)

Step up the number of iterations to 100:

\(\(evenbetter\  = \ 
      Show[realplot, Euler[100], 
\ \ \ \ \ \ \ \ \ \ \ \ 
        PlotLabel -> "100 iterations", 
\ \ \ \ \ \ \ \ \ \ \ \ 
        DisplayFunction -> $DisplayFunction]; \)\)

Step up the number of iterations to 400:

\(\(darngood\  = \ 
      Show[realplot, Euler[400], 
\ \ \ \ \ \ \ \ \ \ \ \ 
        PlotLabel -> "400 iterations", 
\ \ \ \ \ \ \ \ \ \ \ \ 
        DisplayFunction -> $DisplayFunction]; \)\)

Right you are, governor.
The more iterations you use, the better the Euler fake plot you get. 
One more peek:

\(\(Show[
      GraphicsArray[{{rough, better}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {evenbetter, darngood}}]]; 
    \)\)

This is good stuff.

B.1.b) The advantage Euler's method gives you

What is the big advantage you get from Euler's method?

You get a reasonably good plot of a solution even though you might not be able to come up with a formula for the solution.
This is a big, big plus for you.

B.1.c) The fundamental issue of diffeq

In Euler's method, you are given a diffeq
        y'[x]  =  f[x, y[x]] 
       with y[0]  =  b, a given number
You run Euler's and get a nice approximate plot of the solution starting at {0, b}.
How does this illustrate the fundamental issue of diffeq?

In the words of the great American mathematician Peter Lax:
"Every differential equation has a solution and . . . this solution is uniquely determined by initial data [on y[0]]."

B.1.d) The idea behind Euler's faker

Try to explain why Euler's method works.

Here is a differential equation 
         y'[x]  =  Sin[x] y[x]
         with  y[0]  =  2 
and a formula for its solution:

starter = 2; Clear[x, y];
solution=DSolve[{y'[x]==Sin[x] y[x],y[0] == starter},y[x],x];
y[x_] = 
        y[x] /. solution[[1]]

and a plot for 0 ≤ x ≤ 5:

realplot = 
      Plot[y[x], {x, 0, 5}, PlotStyle -> {{Blue, Thickness[0.02]}}, 
        PlotRange -> All, AspectRatio -> (1/GoldenRatio), 
        AxesLabel -> {"x", ""}];

Now look at the Euler faker coming from 8 iterations:

a = 0; b = 5; y[a] = starter; beginner = {a, y[a]}; Clear[x, y, next, jump, iterations, point, k];
jump[iterations_] := ((b - a)/iterations);

next[{x_, y_}, iterations_]:={x, y}+jump[iterations] {1,Sin[x] y};

point[0, iterations_] = beginner;
point[k_, iterations_] := 
        point[k, iterations] = 
          N[next[point[k - 1, iterations], iterations]];
	
Clear[Euler];
Euler[iterations_] := 
        Show[Graphics[{Thickness[0.01], Red, 
              Line[Table[point[k, iterations], {k, 0, iterations}]]}], 
          PlotRange -> All, Axes -> True, 
          AxesOrigin -> beginner, AxesLabel -> {"x", ""}, 
          PlotLabel -> iterations " iterations", 
          DisplayFunction -> Identity];
veryrough = 
        Show[realplot, Euler[8], PlotLabel -> "8 iterations", 
          DisplayFunction -> $DisplayFunction];

Not a very good fake, but this plot exposes what the faker does. 
You can see that the faker actually moves by generating a list of points and moving on straight line segments that connect the points. Here are the points:

\(Clear[fakepoints, fakepointplot]; 

    fakepoints[iterations_] := 
      Table[point[k, iterations], {k, 0, iterations}]; 

    \(fakepointplot[iterations_] := 
      ListPlot[fakepoints[iterations], PlotStyle -> {PointSize[0.03]}, 
        AxesOrigin -> {0, 0}, DisplayFunction -> Identity]; \)\)Show[realplot, fakepointplot[8]];

The faker connects these points with line segments:

Show[realplot, Euler[8], fakepointplot[8], 
      PlotLabel -> "8 iterations"];

When you step up the iteration number, you shorten the jump between successive points:

Show[realplot, Euler[20], fakepointplot[20], 
      PlotLabel -> 20 iterations];

Show[realplot, Euler[50], fakepointplot[50], 
      PlotLabel -> 50 iterations];

Show[realplot, Euler[200], fakepointplot[200], 
      PlotLabel -> 200 iterations];

So far, so good. 

The real question that remains is: 
Why does the fake curve get better and better as you increase the number of iterations?
To answer this, look at the code that generates the fake points:
    next[{x, y}, iterations] = 
       {x, y}  +  jump[iterations] {1,  Sin[x] y};
    point[k, iterations] = 
        N[next[point[k - 1, iterations], iterations]].
This code tells you that once the Euler faker has settled on a point {x, y}, then it moves off that point with the same growth rate as the solution of the differential equation
       y'[x]  =  Sin[x] y[x] 
that passes through that point. So the big differences between the true plot and the fake plot is: 
→ As x advances, the actual solution y[x] updates its growth rate y'[x] instantaneously.  
→ But the fake waits for a while to update the growth rate, updating its growth rate only at each fake point.

When you use a small jump, then you can expect the growth rate of the fake to be nearly the same as the actual growth rate of the function.
To get a small jump on a plotting interval [a, b], you use a high iteration number because 
     jump[iterations] = ((b - a)/iterations).
This is why high iteration numbers usually give better results than low iteration numbers.

B.1.e) The person behind the idea

Say a little bit about the great Swiss mathematician Leonhard Euler, 1707-1783, who originally came up with the idea of faking the plots of solutions of differential equations.

Euler (pronounced "Oiler") was certainly one of the greatest mathematical scientists of all time. In addition to his advanced work in geometry, he wrote the first accurate and complete treatise on calculus. Along the way, he introduced the symbols 
         e, π,  and i  = (Sqrt[-1])
and the abbreviations for Sine, Cosine, Tangent, and the other trig functions.  Even the notation f[x] owes its life to the fertile mind of Euler.
Over and above his mathematics, he had time to take up physics, astronomy, hydrodynamics, optics, electricity, magnetism, ship building, geographical maps, medicine, theology, and oriental languages. Even though he became blind 12 years before his death, nearly half his published works were written during the time he was blind. Not one to neglect the social side of life, Euler spent many years with friends in the famous court of Catherine I in St. Petersburg, Russia.

B.2) Euler's faker and Mathematica's faker, NDSolve

Programmed into Mathematica is an instruction called NDSolve that will accomplish for you everything that Euler's method did for faking the plot of a solution of a differential equation.  And there's an added plus:
You don't have to screw around with setting the number of iterations.
Here's the NDSolve instruction working to get an excellent fake plot for 0 ≤ x ≤ 6 of the differential equation
         y'[x]  =  3.2 Sin[y[x]]  -  0.3 (x^2)     
         with y[0]  =  0.2:

a = 0; b = 4; starter = 0.2; Clear[solution, x, y, fakey];
solution=NDSolve[{y'[x]==3.2 Sin[y[x]] - 0.3 (x^2),y[0] == starter},y[x],{x, a, b}];
fakey[x_] = 
        y[x] /. solution[[1]];

masterfakeplot = 
        Plot[fakey[x], {x, a, b}, 
          PlotStyle -> {{Blue, Thickness[0.01]}}, 
          AxesLabel -> {"x", "y[x]"}, 
          AxesOrigin -> {a, starter}, PlotRange -> All, 
          AspectRatio -> (1/GoldenRatio), 
          PlotLabel -> "Master forgery"];

The NDSolve instruction is just waiting for you to use it.

B.2.a)

Discuss the mathematical principles behind the NDSolve instruction.

The mathematical principles behind the NDSolve instruction are the same as the mathematical principles behind Euler's method.  Except in the case of the NDSolve instruction, extra features typical of professionally written software are included.
You'll hear more about these extra features later.

B.2.b)

Look at this:

a = 0; b = 4; starter = 0.2; Clear[solution, x, y, fakey];

solution=NDSolve[{y'[x]==3.2 Sin[y[x]] - 0.3 (x^2),y[0] == starter},y[x],{x, a, b}];
fakey[x_] = 
        y[x] /. solution[[1]]

What does the output mean?

This output reflects the fact that NDSolve first produces a bunch of points and then strings them together with an interpolating function -- just as Euler's method does. The formula for this interpolating function is not available, but you can plot it:

masterfakeplot = 
      Plot[fakey[x], {x, a, b}, PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y[x]"}, 
        AxesOrigin -> {a, starter}, PlotRange -> All, 
        AspectRatio -> (1/GoldenRatio), 
        PlotLabel -> "Master forgery"];

B.2.c)

But isn't the actual formula always better than a fake?

It all depends on what you want to do.  
If what you want is a plot of the solution, then the fake should be as good as the actual formula.  
If you want to do theoretical analysis of the true solution, then the actual formula is what you want. The only trouble is that the actual formula is available only in special situations; the fake is nearly always available.

B.2.d)

Does this mean that you should give up on Euler's method?

Hell no!
For serious practical situations, you want something better. But as a theoretical tool to use to unlock more secrets of calculus, Euler's method is just the ticket.

B.3) Reading the diffeq through flow plots

B.3.a.i)

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

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

    \(f[t_, y_] = 0.4\ y\ \((1 - y\/b)\); \)\)diffeq = {y'[t] == f[t, y[t]], y[0] == starter}

Here comes a plot of the solution of this differential equation corresponding to given starter data y[0]  =  0.7:

diffeq

\(Clear[y]; 
starter = 0.7; 
\(endtime = 35; \)\)\(y[t_] = 
      y[t] /. \(NDSolve[diffeq, y[t], {t, 0, endtime}]\)[[1
          ]]; 


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

Now look at this:

scale = 6; \(tangentflow = 
      Table[Arrow[{1, f[t, y[t]]}, Tail -> {t, y[t]}, 
          VectorColor -> Blue, ScaleFactor -> scale, 
          HeadSize -> 1.5], {t, 0, endtime, endtime\/20}]; 


    both = Show[tangentflow, solutionplot, PlotRange -> All, 
        Axes -> True, AxesLabel -> {"t", "y[t]"}, 
        AspectRatio -> 1\/2]; \)

Describe what you see and explain why you see it.

Take another look:

Show[both];

What you see is the solution plot being guided by its tangent vectors. This came out this way for the same reason that Euler's method works.
To wit: The vectors plotted are scaled versions of
           {1, f [t, y[t]} with tails at {t, y[t]}.  
Because y[t] is a solution of y'[t]  =  f[t, y[t]], 
these vectors are the same as
           {1, y'[t]} with tail at {t, y[t]}. 
But the slope of 
           {1, y'[t]}  is   y'[t],
which is the instantaneous growth rate of y[t] and that's the reason these are tangent vectors.

B.3.a.ii) Flow plots and corks

Look at this:

scalefactor = 1.1; \(flowplot = 
      Table[Arrow[scalefactor\ {1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.8], {t, 0, endtime, endtime\/15}, {y, 0, 1.4\ b, 
          b\/8}]; 

Show[solutionplot, flowplot]; \)

This shows the solution plotted in part i) going with the flow.
What else can you see from this plot?

Take another look:

Show[solutionplot, flowplot];

Think of it this way:
Each plotted vector is tangent to the solution plot that goes through its tail. Think of the whole t y-plane as fluid flowing in the directions indicated by the tangent vectors plotted above. The plot of a solution is what you get when you drop a cork into the flow and let it float along. 
This way you can use your eyes to see how solutions plot out:

Show[flowplot, Axes -> True];

When you start with 
         y[0]  =  starter between 0 and 10, 
the corresponding solution increases and eventually snuggles onto a straight line. When you start with 
         y[0]  =  starter bigger than 10, 
then the corresponding solution decreases and eventually snuggles onto the same straight line. 
Check it out for a couple of other starting points:

\(Clear[y]; 
diffeq\)

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


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

See these solutions go with the flow:

Show[solutionplots, flowplot];

Math happens.
Play with other starting values.

B.3.b.i) More flow plots

Here's a new differential equation:

\(Clear[diffeq, y, t, f, starter]; 
b = 10; 
\(f[t_, y_] = y - 0.5\ t; 
    \)\)diffeq = {y'[t] == f[t, y[t]], y[0] == starter}

And its flow plot:

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


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


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

Interpret this flow plot.

Take another look:

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

From this flow plot, you can spot at least two types of solutions of the differential equation:

diffeq

Here is a sample of each shown with the flow plot:

\(Clear[y1, y2, y, t]; 
{starter1, starter2} = {0.2, 0.6}; 

    \(endtime = 4.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
          ]]; \(solutionplots = 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; 


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

B.4) Systems of interacting differential equations: 
        The predator-prey model

This  predator-prey  model  was  orginally  proposed  and  studied  by  Volterra  in  1926  in  an  effort  to  explain  the  oscillatory  levels  of  certain  fish  harvests  in  the  Adriatic  Sea.

Here's the idea: 
Two species coexist in a closed environment. One species, the predator, feeds on the other, the prey. There is always plenty of food for the prey, but the predators eat nothing but the prey.
Put 
               pred[t]  =  population  of  predators  at  time  t.
Put  
               prey[t]  =  population  of  prey  at  time  t.

It's reasonable to assume that there are positive constants a and b such that:
               prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
because the abundance of food for the prey allows the birth rate of the prey to be proportional to their current number, and the death rate of prey is proportional to both the current number of prey and the current number of predators.

It also makes some sense to assume that there are positive constants c and d such that:
                pred'[t]  =  -c pred[t]  +  d pred[t] prey[t]
because it's reasonable to assume that the death rate of the predators is likely to be proportional to the current population of predators and that the birth rate of the predators is proportional to both the current number of the predators and the size of the food supply (the prey). This gives you two simultaneous differential equations
                 prey'[t]  =  a prey[t]  -  b prey[t] pred[t], 
                 pred'[t]  =  -c pred[t]  +  d pred[t] prey[t].

B.4.a)

Examine what happens for the sample choice of proportionality constants 
         a  =  0.7,  b  =  0.3,  c  =  0.44,  and  d = 0.08
with the prey starting off with a population of 4.0 units and the predators starting out with a population of 1.2 units. 

Here the units could be thousands or millions so that a fractional population can make sense.

Here is a plot of Mathematica's faker of the prey population as a function of time t for the first 50 time units:

endtime = 50; preystarter = 4.0; predstarter = 1.2; a = 0.7; b = 0.3; c = 0.44; d = 0.08; Clear[pred, prey, t];
approxsolutions=NDSolve[{prey'[t]==a prey[t] - b prey[t] pred[t],pred'[t]==-c pred[t] + d pred[t] prey[t],prey[0] == preystarter,pred[0] == predstarter},{prey[t], pred[t]},{t, 0, endtime}];
pred[t_] = 
        pred[t] /. approxsolutions[[1]]; 
      \(prey[t_] = 
      prey[t] /. approxsolutions[[1]]; 

    
preyplot = 
      Plot[prey[t], {t, 0, endtime}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AxesLabel -> {"t", "prey"}, 
        AspectRatio -> 1\/GoldenRatio]; \)

Hey! A cyclical pattern.
And the predator population as a function of time t:

predatorplot = 
      Plot[pred[t], {t, 0, endtime}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"t", "predators"}, 
        AspectRatio -> (1/GoldenRatio)];

Another cyclical pattern. Here they are together:

\(both = 
      Plot[{prey[t], pred[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.01]}}, 
        AxesLabel -> {"t", ""}, 
        AspectRatio -> 1\/GoldenRatio, 
        DisplayFunction -> Identity]; 


    Show[both, DisplayFunction -> $DisplayFunction]; \)

Look at those populations interacting with each other.

B.4.b)

Here are the plots of both populations shown together again:

Show[both, DisplayFunction -> $DisplayFunction];

Discuss the relationships between the curves.

Cyclic oscillations of both populations are evident as all get out.
Feast your eyes on the vivid relationship between the crests and dips of the predator population and those of the prey population! The cycles are slightly out of sync. When you think about it, they should be slightly out of sync. Look again:

Show[both, DisplayFunction -> $DisplayFunction];

When there are enough prey to sustain strong predator growth, the predators begin to grow so much that they outstrip their food supply. This puts both populations into a decline until the point at which the predators are not a powerful menace; then the prey become numerous enough to support strong predator growth. The cycles go on and on.

B.4.c)

Estimate the length of one cycle in time units. 
Plot one cycle, two cycles, three cycles, and four cycles.

Put in horizontal lines shooting out from the vertical axis reflecting the information that 
         pred[0]  =  0.9 and prey[0]  =  4.

bothwithlines = 
      Plot[{prey[t], pred[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.01]}}, 
        AxesLabel -> {"t", ""}, 
        AspectRatio -> (1/GoldenRatio), 
        Epilog -> {Line[{{0, predstarter}, {endtime, predstarter}}], 
            Line[{{0, preystarter}, {endtime, preystarter}}]}];

A closer look:

Show[bothwithlines, PlotRange -> {{0, 15}, Automatic}];

One cycle is about 12 time units by eyeball estimate.
Use Mathematica to improve on this:

FindRoot[prey[t] == preystarter, {t, 12}]

Check:

FindRoot[pred[t] == predstarter, {t, 12}]

The eyeball was not so bad.
A plot of one cycle:

cycle = 12.2355; one = 
      Plot[{prey[t], pred[t]}, {t, 0, cycle}, 
        PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.005]}}, 
        AxesLabel -> {"t", ""}, 
        PlotRange -> {{0, 4 cycle}, Automatic}, 
        AspectRatio -> (1/GoldenRatio), 
        PlotLabel -> "one cycle"];

Two cycles:

two = 
      Plot[{prey[t], pred[t]}, {t, 0, 2 cycle}, 
        PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.005]}}, 
        AxesLabel -> {"t", ""}, 
        PlotRange -> {{0, 4 cycle}, Automatic}, 
        AspectRatio -> (1/GoldenRatio), 
        PlotLabel -> "two cycles"];

Three cycles:

Plot[{prey[t], pred[t]}, {t, 0, 3 cycle}, 
      PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.005]}}, 
      AxesLabel -> {"t", ""}, 
      PlotRange -> {{0, 4 cycle}, Automatic}, 
      AspectRatio -> (1/GoldenRatio), 
      PlotLabel -> "three cycles"];

And four cycles:

Plot[{prey[t], pred[t]}, {t, 0, 4 cycle}, 
      PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.005]}}, 
      AxesLabel -> {"t", ""}, 
      PlotRange -> {{0, 4 cycle}, Automatic}, 
      AspectRatio -> (1/GoldenRatio), 
      PlotLabel -> "four cycles"];

Nice excursion into mathematical ecology. 
As William E. Boyce and Richard DiPrima point out in their book, Elementary Differential Equations (Wiley, 1977), "Cyclic variations of predator and prey as [predicted above] are often observed in nature." 
In fact E.P. Odum [Fundamentals of Ecology, Saunders, 1953] was able to use the fur catch records of the Hudson Bay Company from 1850 to 1930 to estimate the Canadian lynx and snowshoe hare populations during these years. Odum's curves are not too different in character from the curves you see above.
All mathematical ecologists agree that the predator-prey model is not an end in itself, but is important as a tool in the quest for asking the right questions.

B.4.d)

Why all this fooling around with fake plots?
Wouldn't actual formulas for pred[t] and prey[t] be more convenient?

It all depends on what you want to do.  
If what you want is a plot, then the fake should be as good as the actual formula.  
If you want to do theoretical analysis of the true solution, then the actual formula is what you want. The only trouble is that the actual formula is available only in special situations; the fake is always available.
In the case of the predator-prey model, the formulas so treasured in old-fashioned math classes are not available.
Try it:

Clear[a, b, c, d, pred, prey, t];
DSolve[{prey'[t]==a prey[t] - b prey[t] pred[t],       pred'[t]==-c pred[t] + d pred[t] prey[t],prey[0] == 4,pred[0] == 1},{prey[t], pred[t]},t]

Mathematica  failed to come up with formulas for pred[t] and prey[t]. You can't chalk this up as a fault of Mathematica's, because no scientist has ever succeeded in finding formulas for pred[t] and prey[t]. For the predator-prey model, the fake is the only game in town.

B.4.e)

Where do you go for further reading on predator-prey models?

Predator-prey models have been (and still are) under heavy study by biologists and mathematicians. For further reading, here are some good places to begin:

B.4.f)

What version of Euler's method underlies Mathematica's fake plots above?

To fake a plot of the solution of y'[x]  =  f[x, y[x]], Euler's method uses
        next[{x, y}]  =  {x, y}  +  jump {1,  f[x, y]}
        point[k]  =    N[next[point[k - 1]]]

To fake a plot of the solution of the simultaneous differential equations 
          x'[t]  =  f[t, y[t], x[t]]
          y'[t]  =  g[t, x[t], y[t]], 
Euler's method uses:
 next[{t, x, y}]  = 
         {t, x, y}  +  jump {1,  f[t, x, y], g[t, x, y]}
 point[k]  =  point[k]  =  N[next[point[k - 1]]].
       
The fake points for x[t] as a function of t are fished out from the first two slots of point[k], and the fake points for y[t] as a function of t are fished out from the first and the third slots of point[k].
You could program it yourself, but why bother?