@Carl Love Here is a Maple worksheet with all the computations needed, and the final plots, as I outlined above. The plots match those from the paper. You should put a proper citation for the paper in here, and fill in the descriptions and units of the parameters. The plots could be dressed up a bit, but they clearly show that the computations are accurate.
 

Download Chemoimmunotherapy.mw

@tomleslie s(t) is the input function (or forcing function) of immunotherapy. Typically, it's a square wave, but for the purpose of generating plots (a) and (b), it's constant with value s_infinity. (Note: s0 <> s(0), unless by some extraordinary coincidence.) Likewise for q(t), which is the forcing function of chemotherapy.

@acer Vote up. So, I was wrong: It is possible (to some extent)^[*1] to control the order of multiplicands within a `*` term. I did know that it was possible to use sort to control the order of terms in a `+` expression, but I didn't know that it also worked within the `*`s. Please let me know if I was wrong about anything else.

Would you please show some dismantle examples, as you hinted at, with some exposition? To start with, the distinction between SUM (aka `+`) and POLY seems capricious: 3*a*b*c + a*b*d is a POLY and 3/4*a*b*c + a*b*d is a SUM. Also, I am worried about increasing divergence between the type system and the low-level storage. And by your aside about dismantle above, did you mean that certain data structures in the kernel exist primarily to affect the GUI display at higher levels of prettyprint (or perhaps typesetting)? 

[*1] I say "to some extent" because I still don't think that it's possible to put numeric coefficients any place other than the beginning or to have two expressions with equal `*` terms in different orders.

I'm just restating your Question for those who might find your rambling style confusing:

Can DocumentTools:-Components:-VideoPlayer (or any similar command) be used to stream live video?

Your plots (a) and (b) have their vertical axes labelled D[chemo]. How is D[chemo] determined in terms of the quantities that you've already defined, such as Q(t)?

@Annonymouse Apparently, axis[1]= [mode= log] doesn't work very well when there are bars whose lengths---before taking the logarithm---are so short that they wouldn't appear on a non-log plot. So, we need to take the logs before plotting and then rescale the horizontal axis after. So, here's a better, hand-rolled logarithmic bar chart. Use this after you've assigned the list difference.

difs:= remove(x-> rhs(x)=0, (lhs=abs@rhs)~(difference));
n:= nops(difs);
M:= floor(min((log[10]@rhs)~(difs)));

plot(
   [seq([[0,k], [log[10](rhs(difs[k]))-M, k]], k= 1..n)],
   axis[1]= [tickmarks= [seq(k= typeset(10^nprintf("%d", k+M)), k= 1..-M, 2)], gridlines],
   axis[2]= [tickmarks= [seq(k= typeset(lhs(difs[k])), k= 1..n)]],
   thickness= 9, axes= boxed, view= [default, 0..n+1]
);

Let me know if you need any adjustments.

@Annonymouse As you've no doubt figured out on your own (since apparently you found my Answer useful), axis[2] refers to the vertical axis and axis[1] refers to the horizontal axis. I'm sorry that I didn't catch at first that you were using Statistics:-BarChart, which reverses (IMO) the usual roles of the two axes, so my Answer should've used axis[1].

Still, I'm not satisfied with the results. I'll come up with something that uses regular plot instead.

@Mariusz Iwaniuk Amazing simplification; vote up! How did you guess that useint would help?

My guess is that any "function-like" expression which doesn't contain hypergeom or MeijerG is considered simpler than any that does. By "function-like", I mean to exclude expressions with unevaluated limits, integrals, sums, etc. Those are essentially "meta-functions" to my mind. Because of the way that the arguments of hypergeom and MeijerG are grouped into lists, these functions have an arbitrary number of parameters. Every reduction of the number of parameters could be considered to add a whole "order" of "simplicity". If you look at the definitions of hypergeom and MeijerG on their help pages, you may understand what I mean.

Your intuition is that any proper (or idealized) simplifier S must impose a partial order simpler on its set of possible inputs, at least implicitly, and that S(e) should be minimal wrt simpler. Clearly, this must be true wrt the size option if it's to be meaningful. I think that this partial ordering concept is also true (at least in theory) when simplify is used without an option. The question then becomes What exactly is the partial ordering used? This would be a wonderful addition to the documentation, although it may be the size of a small book.

One consequence of this is idempotency: For any e, S(S(e)) = S(e). I have occasionally (extremely rarely (like twice in 20 years)) seen cases where simplify (with no options) violates this, and these can be considerd bugs.

@nm If f is any symbol, and t is any type, then 'f(t)' is a type that matches functions named f having exactly one argument, with that argument being type t. The difference between 'f(anything)' and specfunc(f) or specfunc(anything, f) is that the latter two will match functions with any number of arguments. Note that the quotes are often needed on 'f(t)' because f is often a procedure name.

This also explains why it didn't complain about your misspelled type.

Christian's Answer seems obfuscated.

@nm You mean extra c, not extra s.

It is interesting that the first 5 of the 6 solutions from dsolve are of the form A_i*(x^(5/2))^(3/5), where the A_i are 5 specific algebraic-number constants. These 5 have no general integration constant. The 6th solution, when solved for y(x), becomes

y(x) = RootOf(_C1*_Z^2*x^(3/2)+x^(5/2)+_Z^5)^3

So, that's essentially the roots of a 5th-degree polynomial.

@syhue I don't understand "ady". Does it work, yes or no?

@Rouben Rostamian  Yes, thanks, that is an improvement.

@Oliveira The uses P= Physics form is safer, but it isn't 100% safe. If you actually execute gds, and it does use P:-diff, the global state will be altered, and it won't be reset when the procedure exits. This global state alteration is of the form shown by Rouben and in your worksheet below. But in neither case will references to diff outside the procedure invoke Physics:-diff rather than ordinary diff.