Mathematica code for the examples of Section 1.13
David McClendon

Compartmental models example 3 (the mixing problem with the salt water tank):

Slope field for the ODE:

phi[t_,y_] :=.6 - .006 y; 
VectorPlot[{1,phi[t,y]},{t,-.5,1000},{y,-.5,200}, VectorPoints -> 35, Axes -> True, VectorScale->{Automatic,Automatic,None}, VectorStyle -> Orange, StreamPoints ->{{0,0}}, StreamScale -> Full, StreamStyle -> {Blue, Thick}]

Euler’s method code:

euler[f_,{t_,t0_,tn_},{y_,y0_},steps_]:=Block[{told=t0,yold=y0,thelist={{t0,y0}},t,y,h},h=N[(tn-t0)/steps];
Do[tnew=told+h;
ynew=yold+h*(f/.{t->told,y->yold});
thelist=Append[thelist,{tnew,ynew}];
told=tnew;
yold=ynew,{steps}];
Return[thelist];]

euler[.6 - .006 y, {t, 0,200}, {y, 0}, 100][[101]]

ListPlot[euler[.6 - .006 y, {t, 0,200}, {y, 0}, 100]]

Solution of the ODE:

g[t_] = 100 - 100 (E^(-.006 t));
Plot[g[t], {t, 0, 1000}, PlotRange -> {0, 150}]

New question:  what is the amount of salt in the tank when t = 200?

New question:  when does the concentration in the tank reach 0.05 kg/L ?

Newtonian mechanics example (falling object):

Slope field for the ODE:

phi[t_,y_] :=9.8 - (3/20)y; 
VectorPlot[{1,phi[t,y]},{t,-.5,80},{y,-.5,100}, VectorPoints -> 30, Axes -> True, VectorScale->{Automatic,Automatic,None}, VectorStyle -> Orange, StreamPoints ->20, StreamScale -> Full, StreamStyle -> {Blue, Thick}]

Euler’s method code:

euler[f_,{t_,t0_,tn_},{y_,y0_},steps_]:=Block[{told=t0,yold=y0,thelist={{t0,y0}},t,y,h},h=N[(tn-t0)/steps];
Do[tnew=told+h;
ynew=yold+h*(f/.{t->told,y->yold});
thelist=Append[thelist,{tnew,ynew}];
told=tnew;
yold=ynew,{steps}];
Return[thelist];]

euler[9.8 - (3/20)y, {t, 0,10}, {y, 0}, 100][[101]]

ListPlot[euler[9.8 - (3/20)y, {t, 0,10}, {y, 0}, 100]]

Solution of the ODE:

Manipulate[Plot[ (196/3)  + C (E^(-3t/20)), {t, 0 , 100}, PlotRange -> {-20,100}], {{C,-65.333},-200,200}]

New question:  suppose the object is dropped from a height of 3000 m.  what is the position of the object at time t?

How long does it take the object to hit the ground?

Heating and cooling:  Example 2

Slope field:

phi[t_,y_] := (1/4)( 70 - 10 Cos[((π t)/12)] - y); 
VectorPlot[{1,phi[t,y]},{t,-.05,24},{y,50,100}, VectorPoints -> 30, Axes -> True, VectorScale->{.02,.02,None}, VectorStyle -> Orange, StreamPoints ->{{0,65}}, StreamScale -> Full, StreamStyle -> {Blue, Thick}]

Euler’s method code:

euler[f_,{t_,t0_,tn_},{y_,y0_},steps_]:=Block[{told=t0,yold=y0,thelist={{t0,y0}},t,y,h},h=N[(tn-t0)/steps];
Do[tnew=told+h;
ynew=yold+h*(f/.{t->told,y->yold});
thelist=Append[thelist,{tnew,ynew}];
told=tnew;
yold=ynew,{steps}];
Return[thelist];]

euler[(1/4)( 70 - 10 Cos[((π t)/12)] - y), {t, 0,10}, {y, 65}, 100][[101]]

ListPlot[euler[(1/4)( 70 - 10 Cos[((π t)/12)] - y), {t, 0,24}, {y, 65}, 1000], PlotRange -> {50,100}]

Solution of the IVP:

M[t_] = 70 - 10 Cos[((π t)/12)];
K = 1/4;
f[t_] = ∫(E^(K t)) K M[t]\[DifferentialD]t

Solve[65 == (E^(-0/4))(f[0] + C), C]

T[t_] = (E^(-t/4))(f[t] + -((5 (-9+(π^2)))/(9+(π^2))))
Plot[{M[t], T[t]}, {t, 0, 24}, PlotRange -> {50,90}, PlotStyle -> {{Thick, Red}, {Thick, Blue, Dashed}}]

New question:  what is the temperature inside the house at 11 AM?

New question:  what is the maximum temperature inside the room?  At what time is the room warmest?