To go backwards, switch the from and to and add by -1: Change the second for-loop header to

for i from n-1 to 1 by -1 do ...

Note that your n-2 was also incorrect; it should be n-1. The order of the prepositions for, from, to, and by doesn't matter.

Your ODEs can both be factored, as is shown by you getting multiple solutions when solving for the (highest order) derivative. Perhaps in those cases odeadvisor only analyzes the simpler factor. If there are multiple factors, you could map odeadvisor over them.

This is a start for how to get plots (a) and (b):

Given their similarity, it's no surprise that these two plots contain the same information in different formats. The relationship that connects the two plots is D__chemo = q__infinity * t__cure (total dose = constant administration rate * time administered). In other words, the vertical axis on both plots is the product of the horizontal axes of the individual plots. And t__cure is found by solving N(t) = 1 for t (or perhaps N(t) = 0.5 is a better choice). In other words, the patient is cured when the number of live cancer cells is less than 1.

Outline of steps:

1. Replace s(t) and q(t) with s__infinity and q__infinity in the ODE system.
2. Make these formal parameters to the dsolve command by using option parameters= [s__infinity, q__infinity].
3. Use the events option of dsolve to set a "trigger" for N(t) = 1 or N(t) = 0.5.
4. Call dsolve to return the solution procedure, which I'll call dsol. This only needs to be done once.
5. Write a procedure (I'll call it t_cure) that accepts a numeric value of q__infinity as its parameter and returns the corresponding value of t__cure by passing q__infinity to dsol as a parameter and integrating until the event "trigger" is "fired" and reading the value of t. 
6. Set a value of s__infinity via dsol(parameters= [s__infinity= ...]).
7. Use procedure t_cure to make a semilogplot of t__cure vs. q__infinity over some reasonable interval based on the plot (a) that you want to end up with. The purpose of this semilogplot is just to collect numeric data; there's no need to look at it if your goal is just to get plots (a) and (b).
8. Extract the two-column data matrix from the plot via op([1,1], ...). The first column is values of q__infinity, and the second is t__cure.
9. Multipy the two columns elementwise. These are the values of D__chemo.
10. Go back to step 6. Repeat with different value of s__infinity. 
11. Using regular plot with mode= log on both axes, plot D__chemo vs. q__infinity; this is plot (a). Do likewise for D__chemo vs. t__cure; this is plot (b). 

Add the option axis[2]= [mode= log] to the plotting command.

All that you need to do is to change specfun to specfunc.

Another type that you might find useful is typefunc(mathfunc).

You'd use subsindets if you want to take the items found by indets, modify them, and then substitute the modifications for the originals in the original expression.  It's very powerful.

You are trying to merge a static plot and an animation. A convenient way to do that is the background option to animate. Like this:

plots:-animate(
   plots:-arrow,
   [eval([p, ptan], s= A)[], width= 0.3, length= 4],
   A= 0..3, digits= 10,
   background= plots:-spacecurve(p, s= 0..16*Pi, numpoints= 1000)
);

The session transcripts that you show do not illustrate the phenomenon that you claim. In your first example, the error message about Physics:-diff is coming from your procedure gds which should be using Physics:-diff. In your second example, the unevaluated Determinant is external to procedure dst.

There is no bug illustrated here, nor even undesirable behavior, regardless of whether you use sec or seq. Rouben's example, however, does show a different form of undesirable behavior---a change in the global state.

Here's the worksheet:
 

Download RLC.mw

You made a small mistake while solving. When you divide through by 10, the last term becomes 10*x. You have 50*x.

Here's how to check your work with Maple:
 

Download RLC.mw

This will work for both sums and products, and many other types as well:

op(0, expr)(op(2.., expr))

The syntax M[i][j] is supported (although it can be extremely inefficient), and your assignment statements are actually being executed. So something(s) is/are being assigned the value 1. The questions are Why doesn't this produce your expected result? And what happened to those somethings? I submit the following two examples as aids towards you understanding that, knowing that you are a sophisticated user. I do not suggest that you or anyone else do this in actual polished code.

doit:= proc()  
local i, j, Row, M:= Matrix((3,3));                                                         
   for i to 3 do
      Row[i]:= M[i];  
      for j to 3 do  
         Row[i][j]:= 2^i*3^j                                     
      od;
      M[i]:= Row[i]                             
   od;              
   M 
end proc
:                                           
doit();  

doit:= proc()  
local 
   i, j, 
   M:= Matrix((3,3), order= C_order), 
   ByRows:= Vector(3, k-> ArrayTools:-Alias(M, 3*(k-1), [3]))
;                                                         
   for i to 3 do  
      for j to 3 do  
         ByRows[i][j]:= 2^i*3^j                                     
      od                             
   od;              
   M  
end proc
:                                           
doit(); 

Both examples produce the same result---a Matrix M equivalent to Matrix((3,3), (i,j)-> 2^i*3^j). (I chose this Godel-numbering formula because you can immediately reconstruct the i and j from the entry alone.) The syntax M[i] extracts the ith row as an anonymous Vector, and M[i][j]:= ... assigns to the jth element of that Vector. What happens to that anonymous Vector? Unless you assign it back to M[i] (as in the first example), it gets garbage collected. The second example uses Alias to rename M's row vectors so that assignment to them (through the name ByRows) is equivalent to assignment to M itself rather than to some anonymous object. So, this second example allows you to use syntax very similar to your original. But, once again, I do not recommend that you actually do this because there's no good reason to since there are several much better ways.

I'm sure that by now you know about unapply, as in unapply(f,x). However, it doesn't work well in this double-arrow situation. Cases that are too complicated for unapply can usually be handled by a special usage of subs:

g:= subs(_f= subs(_f= f, x-> _f), (a,b)-> _f)

Note that unlike Kitonum's Answer, this does produce a true double-arrow procedure:

(a,b)-> x-> a*x+b

This usage of subs is AFAIK undocumented, and it cannot be fully understood through an understanding of ordinary subs. However, it has extensive and long-standing use in library code, such as in the code of unapply itself. For me, it's a fundamental tool of meta-programming, and it's often the only way to accomplish a goal, such as producing a procedure that can be compiled to external code or used with evalhf or Threads. The procedure g produced by Tom's method cannot be used like that.

I just showed you yesterday how to generate 2000-bit primes! Did you see that Answer?

There are no special commands needed to work with large integers in Maple. You just enter the computation normally, for example: 2^2000. It can easily handle numbers even much larger than that.

@Christian Wolinski 

Aha! I managed to deconstruct subsindets. The basic built-in subsindets (by which I mean without any options in square brackets) is equivalent to this Maple procedure: 

SubsIndets:= proc(e, t::type, T)
local J:= sort([indets(e, t)[]], (x,y)-> length(x)<length(y)), r:= e; # [*1]
   to nops(J) do (J,r):= subs(J[1]= T(J[1], _rest), [J[2..], r])[] od;
   r
end proc:

So, applying it to your examples:

e:= (((Q(a)^3)^(5/4))^(15/7))^(6/8);
T:= x-> `if`(op(1,x)::specfunc(Q), Q(x), x);
SubsIndets(e, anything^rational, T);
SubsIndets(e, specfunc(Q)^rational, Q);

So, considering the simplicity of the recursive process in the above algorithm, this may be the intended behavior. If so, then the take-away is that the user should refine the type in the second argument rather than filtering to some subset of that type in the transformation procedure. If you have trouble refining the type, you can always resort to something like this:

And(t, satisfies(proc(x) ... end proc))

[*1] As of Maple 2017 (I believe), that sort command can be shortened to

sort([indets(e,t)[]], key= length)

The commands nextprime or prevprime can do it very quickly. Here's three examples: the smallest 2000-bit prime, the largest, and a random one:

p:= nextprime(2^1999);
q:= prevprime(2^2000);
randomize(): r:= nextprime(rand(p-1..q-1)());