This is Tom's summary of section 1.8 of our textbook, Differential Equations-A Modeling Perspective, by Borrelli and Coleman. You can grab a copy of the actual Maple worksheet by right-clicking the link and selecting "Save Link As" from the ensuing menu.

Load the typical packages:
> map(with, {DEtools, plots, plottools}):

Warning, new definition for translate

By now, this solution should not surprise you:
> x_ans:=dsolve(xeqn, {x(t)});

Now, substitute that into the ODE for y(t). The form of the answer above is ready made for substitution, since I enclosed x(t) in set braces when I solved for x(t).
> new_y:=subs(x_ans, yeqn);

That too is linear and the integrating factor is exp(k[2]*t). The right side integral is the product of two exponentials (not too bad, as integrating C*exp([k2 - k1]*t) is straightforward), but we'll let Maple solve it.
> y_ans:=dsolve(new_y, {y(t)});

This cleans it up a bit (I should really use rhs(y_ans), but this worked without it):
> y_ans:=factor(y_ans);

We're going to want to plot these solutions as expressions in t, and plug in various values for A, k[1] and k[2]. I'll just handle this one case at a time, doing the substitution, then making it into a function (if needed). The first examples use these values of the parameters:
> vals:= {A=1, k[1]=0.6931, k[2]=0.0231};

Insert these values (into the rhs of the answers):
> x1:=subs(vals, rhs(x_ans) );
y1:=subs(vals, rhs(y_ans) );

For plots, we can use expressions. The GI tract is blue and the Bloodstream is red.
> plot([x1,y1], t = 0..6, color = [blue,red],
labels = ["hours","Units"] );

To get several values for k[2], we make the other substitutions, and turn k[2] into a plain k (which is easier to work with):
> y2:=subs({A = 1, k[1] = 0.6931, k[2] = k}, rhs(y_ans) );

Make that into a function of k:
> y2 := unapply(y2, k);

And plot it for various values of k.
> plot([y2(.231), y2(.0731), y2(.0231), y2(.00731), y2(.00231)],
t = 0.. 20, color=red);

The text describes a "square wave" function sqw(t, d, p) (in appendix B). The input variable (time) is t. d and p are parameters which describe the duty (on time percentage) and period (amount of time, t, before the graph repeats). You can view sqw(t,d,p) as a function which is ON (equal to one) for the first d percent of a period, and OFF (equal to zero) for the rest of a period. The periods start at integer multiples of p, and last for p time units (0 to p, p to 2p, 2p to 3p etc.). For a cold pill which is on for the first 1/2 hour of each 6 hour time period, then off for the last 5.5 hours, we can take p = 6 hours, and d = 1/12 (of a period) expressed as a percent, which is 100*1 / 12 = 25 / 3 percent. Here is a definition of the general sqw function, which you can copy and paste and reuse in the future. You should execute this command (once) before attempting to re-use it.
> sqw := (t,d,p) -> piecewise(p*(t/p-floor(t/p)) <= d*p/100, 1, 0);

Here is the texts I(t) = on for 0.5 hours at the start of each 6 hour period.
> plot(sqw(t,25/3,6), t=0..25, thickness=2);

I prefer to think in terms of ONTIME and PERIOD for this assignment, so here is a shortcut to create a square wave function which is on for the first "on" hours of a period of length "period".
> onfor:=(on,period)->sqw(t, 100*on/period, period);

The same plot as above:
> plot(onfor(0.5,6), t = 0..25, thickness=2);

To deliver 6 units of antihistamine (at a steady rate) over 0.5 hours, we deliver antihistamine at 12 units per hour for 1 / 2 hours, or 12*onfor(0.5, 6) as a "rate in" function. Like the text, I'll deligate the work of solving this on-off ODE to Maple's numerical plotter. Here is the system (again, separated into two equations, with no initial conditions, since they are entered separately for DEplot).
> xeqn:=diff(x(t), t) = 12*onfor(0.5, 6) - 0.6931*x(t):
yeqn:=diff(y(t), t) = 0.6931*x(t) - 0.0231*y(t):

The GI plot, bloodstream plot (i.e. the system with a scene showing y versus t), and then display them together. Note the "resetting" of stepsize to a small value.
> giplot:=DEplot({xeqn}, x(t), t = 0..45, [ [x(0)=0] ],
arrows = none, linecolor = blue, stepsize = 0.025):

> bsplot:=DEplot({xeqn, yeqn}, {x(t), y(t)}, t = 0..45,
[ [x(0) = 0, y(0) = 0] ], arrows = none, linecolor = red,
stepsize = 0.025, scene = [t,y] ):

> display( [giplot, bsplot], labels = ["hours", "Units"] );

Here is a fancy way to produce a DEplot for some value of k[2] (say k) with k[1] = 0.6931.
> makeit:=k->DEplot({xeqn, diff(y(t), t) = 0.6931*x(t) - k*y(t)},
{x(t), y(t)}, t = 0..120,
[ [x(0) = 0, y(0) = 0] ],
arrows = none, linecolor = red,
stepsize = 0.1, scene = [t,y] ):

Here are definitions of the "zone lines":
> low:=plot(5, t = 0..120, color = black):
high:=plot(40, t = 0..120, color = black):

And a plot of several bloodstream graphs along with the therapeutic zone lines (this takes quite a while):
> display( [low, high, makeit( 0.0462), makeit(.0231), makeit(.0077) ],
labels = ["hours", "Units"] );

Return to Differential Equations Materials.

Return to Maple Materials Page.