DE.01 Transition from Calculus to DiffEq:
The Exponential Differential Equation
       y'[t]+r y[t]=f[t] 
 BASICS

\(Needs["Graphics`Colors`"]; 

    Graphics`Colors`GosiaGreen = RGBColor[0, \ 0.392187, \ 0]; 

    Needs["Graphics`FilledPlot`"]; 
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

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] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/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] == 2Vector[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[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[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[]; \)\)

VectorField and Phase Lines

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, #1 ≥ a && #1 ≤ b&]; 
        zero1 = Sort[zero ⋃ {a, b}]; len = Length[zero]; 
        len1 = Length[zero1]; 
        signs = Table[
            Sign[fcn /. 
                var -> 
                  (1/2) (
                      zero1[[k + 1]] + 
                        zero1[[k]])], {
              k, 1, len1 - 1}]; 
        stick = {Thickness[0.007], Line[{{pos, a}, {pos, b}}]}; 
        dots = (Point[{pos, #1}]&)/@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 slen)/4)]; 
        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, #1 ≥ a && #1 ≤ b&]; 
        zero1 = Sort[zero ⋃ {a, b}]; len = Length[zero]; 
        len1 = Length[zero1]; 
        signs = Table[
            Sign[fcn /. 
                var -> 
                  (1/2) (
                      zero1[[k + 1]] + 
                        zero1[[k]])], {
              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]], #1}]&)/@
            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), (slen/2)]; 
        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, #1]&)/@Options[picture]; 
          opt3 = opt3 //. {
                (AxesLabel -> _) -> AxesLabel -> None, 
                (PlotRange -> _) -> 
                  PlotRange -> Automatic, 
                (PlotRegion -> _) -> 
                  PlotRegion -> Automatic, 
                (Ticks -> _) -> Ticks -> Automatic}; 
          picture = List@@picture; 
          picture[[-1]] = opt3; 
          picture = 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[];

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, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

B.1) The unforced exponential differential equation, the most important differential equation of all: 
                   y'[t] + r y[t] = 0.
This is the boss diffeq.
                     y[t] =  (E^(-r t)) starter
solves the unforced exponential diffeq
          y'[t]+r y[t]=0 with y[0] = starter.

B.1.a)    y[t]  =  a (E^(-r t))solves  y'[t]+r y[t] = 0

Given a number r, to solve the exponential differential equation 
       y'[t]  +  r y[t]  =  0, 
you want a function y[t] with the property that 
       y'[t]  +  r y[t]  =  0  for all t's.
How do you get it?

You do educated guesswork.
The diffeq  is
       y'[t] + r y[t]=0 for all t's.
You go with numbers K and z and put 
       guessy[t] = K (E^(z t))
Because the diffeq is
       y'[t]  +  r y[t]  =  0, 
you want to select z so that 
       guessy'[t] + r guessy[t]=0.
This is the same as
       K z (E^(z t))  +  r K (E^(z t )) =  0.
 This is the same as
       z = -r 
So 
       guessy[t] = K (E^(-r  t)).
Check it out:

\(Clear[guessy, K, t, r]; 
guessy[t_] = K\ E\^\(\(-r\)\ t\)\)

Check whether
       guessy'[t] + r guessy[t] = 0:

guessy'[t]+r guessy[t]

Good.
This tells you and everyone else that
       y[t] = a (E^(-r  t))
solves the exponential diffeq
       y'[t]+ r y[t] = 0.

You can get Mathematica to spit this out directly:

\(Clear[y, t, r]; 
DSolve[\(y'\)[t] + r\ y[t] == 0, y[t], t]\)

This means C[1] can be any constant. You can reconcile the solution 
        y[t] = K (E^( -r  t)) 
above with the Mathematica solution 
        y[t] = C[1] (E^( -r  t))
by taking K =  C[1].

B.1.b)

Solve the diffeq
             y'[t]+0.6 y[t] =0 with y[0] = 12.8 . 
Give a nice plot of the solution.

You already know
          y[t]   =  K (E^( -0.6 t)) 
 solves
          y'[t]+ 0.6 y[t] = 0.
         
But when you look at y[0], you get:
         y[0] = K (E^0) = K.
This tells you that to solve
          y'[t]+ 0.6 y[t] = 0 with y[0] = 12.8, 
all you have to do is to write down
           y[t]  = 12.8  (E^( -0.6 t)).
 Check it out:

\(r = 0.6; 
starter = 12.8; 
Clear[y, t]; 

    y[t_] = starter\ E\^\(\(-r\)\ t\)\)

See this solution:

\(endtime = 20; 

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

It falls off to 0 mighty fast. 
That's what lots of folks like to call "exponential decay."

B.1.c.i)

Given the diffeq
             y'[t]+ r y[t] = 0 with y[0] = 7.6, 
 how does the sign of r tell you whether the plot of the solution goes up or goes down to 0?

If r < 0, solutions go up. If r > 0, solutions go down.

Reasons:
If r < 0, then you get 
       y'[0] + r y[0]=0
or      
       y'[0]=-r y[0]
       y'[0]=-7.6 r 
Because r < 0 tells you that y[t] goes up as t advances from 0.
Then when you look at
       y'[t]=-r y[t],
you see that as y[t], goes up so does y'[t].
The upshot:
       y[t] grows at an increasingly rapid rate.
This increasingly rapid growth is called exponential growth.
See the solution of y'[t] + r y[t] = 0 with y[0] = 7.6
with r = -0.8:

\(starter = 7.6; 
r = \(-0.8\); 

Clear[y, t]; 

    y[t_] = starter\ E\^\(\(-r\)\ t\)\)\(endtime = 20; 

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

Look at that bastard grow. That's big time exponential growth.

On the other hand, if r > 0, 
the solution is y[t] = 7.6 (E^(-r t)) which is guaranteed to decay to 0 pretty durn fast:
See the solution of y'[t] + r y[t] = 0 with y[0] = 7.6 and
with r = 0.8:

\(starter = 7.6; 
r = 0.8; 

Clear[y, t]; 

    y[t_] = starter\ E\^\(\(-r\)\ t\)\)\(endtime = 20; 

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

Decaying to 0. Typical exponential decay.

B.2)   y[t] =  (E^(-r t)) starter +   (E^(-r t))(∫_0^t) (E^(r t)) f[t] dt
solves the forced exponential diffeq
          y'[t]+r y[t]=f[t] with y[0] = starter.

B.2.a)

Old hands at diffeq know that you can get a formula for the solution of the forced exponential diffeq
       y'[t]  +  r y[t]  =  f[t], with y[0] = starter
merely by calculating an integral. 
The formula for the solution is
       y[t] = (E^(-r t)) starter +   (E^(-r t))(∫_0^t) (E^(r t)) f[t] dt . 
            
Here is the formula at work on 
       y'[t]+2.6 y[t]=3.4  (E^( -0.3 t)), with y[0] = 5.1:

\(Clear[yformula, f, t]; 
f[t_] = 3.4\ E\^\(\(-0.3\)\ t\); 
r = 2.6; 

    starter = 5.1; 

    yformula[t_] = 
      E\^\(\(-r\)\ t\)\ starter + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\)\)

See it:

\(endtime = 10; 

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

Here is the formula at work on 
       y'[t]-1.2 y[t]=0.8  t, with y[0] = 0:

\(Clear[yformula, f, t]; 
f[t_] = 0.8\ t; 
r = \(-1.2\); starter = 0; 

    yformula[t_] = 
      E\^\(\(-r\)\ t\)\ starter + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\)\)

See it:

\(endtime = 8; 

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

Here is the formula at work on 
       y'[t]+ 0.5 y[t]=3.7 Sin[2 t], with y[0] = -3.0:

\(Clear[yformula, f, t]; 
f[t_] = 3.7\ Sin[2\ t]\)\(r = 0.5; 
starter = \(-3.0\); 

    yformula[t_] = 
      E\^\(\(-r\)\ t\)\ starter + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\)\)

See it:

\(endtime = 20; 

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

Explain why the formula
       yformula[t]  =  
        (E^(-r t)) starter +  
          (E^(-r t)) 
            (∫_0^t) (E^(r t)) f[t] 
                dt
is guaranteed to work (if the integral can be done)

Start with
       y'[t]  +  r y[t]  =  f[t], with y[0] = starter.
Multiply both sides by  (E^(r t)) to get:
       (E^(r t)) y'[t]+ (E^(r t)) r y[t]=(E^(r t)) f[t].

And you note (thanks to the product rule of derivatives) that the new equation is same as
       D[(E^(r t))  y[t],  t]  =  (E^(r t)) f[t].
 Check:

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

    ∂\_t\((E\^\(r\ t\)\ y[t])\) == E\^\(r\ t\)\ f[t]\)

This tells you that the original diffeq 
       y'[t]  +  r y[t]  =  f[t]
 is the same as 
       D[(E^(r t))  y[t],  t]  =  (E^(r t)) f[t].
 Next you integrate both sides  of
       D[(E^(r t))  y[t],  t]  =  (E^(r t)) f[t]
from 0 to t to get
       (E^(r t)) y[t] -  (E^0) y[0] 
        = 
        (∫_0^t) D[((E^(r t)) y[t]), t] dt = 
          (∫_0^t) (E^(r t)) f[t] dt.
This gives         
       (E^(r t)) y[t] - y[0] = 
        (∫_0^t) (E^(r t)) f[t] dt.
 So
       y[t]  -  (E^(-r t)) y[0] = 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] dt

This gives the exact formula for y[t] in terms or f[t], r and y[0]:
       y[t] =  
        (E^(-r  t)) y[0] +  
          (E^(-r  t)) 
            (∫_0^( t))(E^(r  t)) f[t] 
                dt.

B.2) The downside

What's the downside of trying to integrate out?

If you're working by hand with pencil and paper, the integrals can drive you nuts.
It can even drive mighty Mathematica nuts

Here is the formula at work on 
       y'[t]+ 0.5 y[t]=3.7 Sin[2 t], with y[0] = -3.0:

\(Clear[yformula, f, t]; 
f[t_] = 3.7\ E\^Sin[t]\)\(r = 0.5; 
starter = \(-3.0\); 

    yformula[t_] = 
      E\^\(\(-r\)\ t\)\ starter + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\)\)

Mathematica could not calculate (∫_0^t)(E^(0.5` t + Sin[t]))dt.
Neither can anything or anybody else.
If you can calculate the integral, the method is as easy as apple pie. 
But if you can't do the integrals, the method is dead in the water.

Moral: Integrating out works best on contrived classroom problems.

B.3) Steady state for y'[x]+r y[x]=f[t] with r > 0 : Any solution eventually settles into the same steady state. If r < 0, then all bets are off

B.3.a.i)

Here are three plots of solutions a random forced exponential diffeq 
       y'[t]+r y[t]=4 Sin[2.5 t],  ( with r > 0)
with random starting values on y[0]:

\(r = Random[Real, {0.5, 1.0}]; 
f[t_] = 4\ Sin[2.5\ t]; 
endtime = 20; 
    

Clear[y, y1, y2, y3, t]; 
starter1 = Random[Real, {5, 10}]; 
    starter2 = Random[Real, {\(-2\), 2}]; 
    starter3 = Random[Real, {\(-10\), \(-5\)}]; 
    y1[t_] = E\^\(\(-r\)\ t\)\ starter1 + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\); 
    y2[t_] = E\^\(\(-r\)\ t\)\ starter2 + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\); 
    y3[t_] = E\^\(\(-r\)\ t\)\ starter3 + 
        E\^\(\(-r\)\ t\)\ 
          \(∫\_0\%t\( E\^\(r\ t\)\ f[t]\) \[DifferentialD]t\); 
    plots1 = Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.018], Blue}, {Thickness[0.014], Red}, {
              Thickness[0.01], CadmiumOrange}}, PlotRange -> All, 
        AspectRatio -> Automatic, AxesLabel -> {"t", ""}, 
        PlotLabel -> "y'[t] + r y[t] = f[t] with r > 0"]; \)

Describe what you see and explain why you see it.

Take a look at another one - this time it will be
       y'[t] + r y[t]=6 (E^(-0.4 t))  +  p t  
with random r > 0 random choice of p and random starting values on y[0]:

r = Random[Real, {0.5, 1.0}];
p = Random[Real, {-0.5, 0.5}];f[t_] = 6 (E^(-0.4 t)) + p t;
endtime = 20;

Clear[y, y1, y2, y3, t];
y'[t]+r y[t]==f[t]starter1 = Random[Real, {5, 10}]; starter2 = Random[Real, {-2, 2}]; 
    starter3 = Random[Real, {-10, -5}]; 
    y1[t_] = (E^(-r t)) starter1 + 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] \[DifferentialD]t; 
    y2[t_] = (E^(-r t)) starter2 + 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] \[DifferentialD]t; 
    y3[t_] = (E^(-r t)) starter3 + 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] \[DifferentialD]t; 
    plots2 = Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.018], Blue}, {Thickness[0.014], Red}, {
              Thickness[0.01], CadmiumOrange}}, PlotRange -> All, 
        AspectRatio -> Automatic, AxesLabel -> {"t", ""}];

Every time the three solutions begin their trip in totally different ways, but when t is large, you need a scorecard to tell them apart.  
In other words, they all settle into the same steady state behavior- regardless of the starting values on y[0]. The long term behaivior of the solutions is totally insensitive to the starting value on y[0].

To explain why this happens, look at the formula
       y[t] = 
        (E^(-r t)) starter + 
          (E^(-r t)) 
            (∫_0^t) (E^(r t)) f[t] 
                dt

If r > 0, then as t gets large, you are guaranteed that (E^(-r t)) starter -> 0.
The upshot: In the long run,  (E^(-r t)) 
        (∫_0^t) (E^(r t)) f[t] dt
That's what you're seeing here:

\(Show[plots1]; 
Show[plots2]; \)

B.3.a.ii)

Does this mean that when you have a a forced exponential diffeq 
       y'[t]+r y[t]=f[t], with r > 0 
 then the long term steady state behavior of solutions  
-> depends only on the behavior of f[t] 
 and
-> does not depend not on the starting value for y[0]?

B.3.a.ii)

What happens when r < 0?

Try it and see:

r = -0.1;
f[t_] = 3.8 Sin[2 t];
endtime = 20;

Clear[y, y1, y2, y3, t];
y'[t]+r y[t]==f[t]starter1 = Random[Real, {5, 10}]; starter2 = Random[Real, {-2, 2}]; 
    starter3 = Random[Real, {-10, -5}]; 
    y1[t_] = (E^(-r t)) starter1 + 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] \[DifferentialD]t; 
    y2[t_] = (E^(-r t)) starter2 + 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] \[DifferentialD]t; 
    y3[t_] = (E^(-r t)) starter3 + 
        (E^(-r t)) 
          (∫_0^t) (E^(r t)) f[t] \[DifferentialD]t; 
    plots = Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.018], Blue}, {Thickness[0.014], Red}, {
              Thickness[0.01], CadmiumOrange}}, PlotRange -> All, 
        AspectRatio -> (1/GoldenRatio), 
        AxesLabel -> {"t", ""}, 
        PlotLabel -> "y'[t] + r y[t] = f[t] with r < 0"];

If r < 0, the solutions don't merge. 
To see why, look at the formula
       y[t] = 
        (E^(-r t)) starter + 
          (E^(-r t)) 
            (∫_0^t) (E^(r t)) f[t] 
                dt.
When r < 0, the term (E^(-r t)) starter does not go away as t gets big.
In fact (E^(-r t)) starter gets bigger and bigger as t gets large. 
That's why different starter values on y[0] yield increasingly different plots as t gets large.