I don't know if this is what you're looking for, but you could figure the prime moduli over which the polynomial factors, like this:

for k to 20 do
   p:= ithprime(k);
   F:= Factor(n^2+n+39) mod p;
   if F::{`*`, `^`} then
      print('p'=p, F)
   fi
od:

1. node_data[1.., 1] extracts the first column from Matrix node_data. It could also be written node_data[.., 1].

2. If A and B are Matrices or Vectors with the same number of rows, then <A | B> is the Matrix formed by joining them horizontally.

3. If S is a sequence, then < S > is the corresponding column vector.

4. add(ad_mat[i, j] * wt_mat[i, j], j = 1.. num_characters) multiplies the ith rows of the two Matrices elementwise and then adds those products. It could also be written add(ad_mat[i, ..] *~ wt_mat[i, ..])

5. seq(..., i = 1. . num_characters) does the above computation for every row.

6.  The line BetweenessCentrality_sorted := FlipDimension( x[2])))>, 1) has unbalanced brackets. It isn't syntactically correct, so it does nothing.

I suppose that you want to solve for {t, x, y, z}? The system is quadratic, not linear, in those variables, but we can proceed anyway. The system cannot be completely solved for those variables, but you can get a conditional solution with eliminate.

eliminate({eq||(1..4)}, {t,x,y,z});

This returns a list of two sets. The first set is equations solved for those 4 variables. When I ran it, I got 

{t = 1/2*(a-b)/z, x = -(j+m)*z/(a-b), y = -1/2*(a^2-b^2)/(j+m)/z, z = z}

which shows z as a free variable. (You may get another variable as the free variable.) The second set is expressions in the other variables ({a, b, j, m} in this case) that must be equal to 0 for the first set to be valid. In this case, I got 

(a-b)*(a+b)*(a^2-b^2+j^2-m^2)

In other words, in order for there to be a solution for {t,x,y,z}, either a = b or a = -b or a^2 - b^2 = m^2 - j^2.

The solve command will not return these conditional (or singular) solutions. 

If you're talking about a boundary value problem (BVP) in ordinary differential equations (ODEs), then variable grid sizing is handled automatically by dsolve. It can be controlled with the options initmesh and maxmesh.

If you're talking about partial differential equations (PDEs), then these are called multigrid methods. These are far beyond what's currently implemented in Maple for the numeric solution of PDEs. You'd need to program your own. It'd be an extensive project, but ultimately doable.

Yes, using vectors is almost necessary. I can think of a convoluted way without vectors, but it's just an academic exercise and definitely not worth the trouble.

If you know that you're going to be exporting to Excel, you should probably be using Vectors from the very start, with datatype= hfloat. This would probably be more computationally efficient than lists. Note that you don't need to know the final number of elements in a Vector while you're building it, and that it's far, far more efficient to build a variable-length Vector than a variable-length list.

If you're using Maple 2019, there's some new syntax that simplifies building variable-length Vectors. Go to the help page ?assignment and read the section "Special Semantics of the ,= Assignment Operator". 

Use the setting

interface(prettyprint= 3):

This is usually set as the default. Somehow yours got changed. (Oddly, this is at least the second time recently that I've seen a user report this.) You can make it the default by using the menu commands Tools => Options => Display => Output Display => (2-D Math Notation or Typeset Notation) => Apply Globally.

The topic is called cluster analysis. You'll find much information on the Web about it. I wrote some Maple code a few years ago to do it. Someone who is better at searching MaplePrimes than I am should be able to find it.

I see nothing wrong in your worksheet. Radical algebraic expressions can appear in numerous equivalent forms, for example 1/sqrt(2) = sqrt(2)/2. Your collect command does nothing in this case: The result of your asympt is a series structure that's automatically collected into powers of k, so collect isn't needed.

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.