You can use PolynomialTools:-Hurwitz to find the conditions for stable roots, i.e for roots with negative real parts.

Download Hurwitz.mw

I also would like an easy way to do this. I think evala(Reduce()) and evala(Normal()) want to have the the answer in the same algebraic field extension so don't do this. It can be done somewhat awkwardly through Minpoly; perhaps you had another objection to this. 

Download RootOf.mw

Download Matrix.mw

See here for simplifying the entry (and also the coversion to the matrix).

The index=real[4] in the RootOf looks wrong. [Edit: I submitted an SCR].

Download solve.mw

Download remove_undefined.mw

If I assign some values such as h_n:=[5,4], etc (for np=2) then it works.

points.mw

or a variation more suitable to a large number of points.

points2.mw

Your variable s__v_0 actually had an unseen multiplication in. I replaced it with just sv0. Some subscripts were made differently in different places, for example v[1] or v__1; I changed them all to the double-underline form. In your differential equation mu was strange (from a font perhaps instead of typing mu and using escape), and different from the mu in the parameters list. You had two alphas in the parameter list.

You now need to resolve the fact that yor equations have Q(t) but there is no differential equation for Q(t). You have some parameters in your parameter list that are not in your equations - I would have removed them but you give them values later. So next step is to resolve these inconsistencies.

Covid19_Simulation2.mw

Alias perhaps comes closest to "doing it automatically".

Download alias.mw

I just removed all entries with bonds with H from the table, then made the matrix, rather than deleting rows and columns afterwards; I hope this is equivalent.

graph:=proc(t::table)
  local t2,n,A;
  t2:=remove['flatten'](x->x[1][-1]="H",t);          # remove all entries from table with string ending in "H"
  n:=max(ListTools:-Flatten([indices(t2)]));                   # find Matrix size
  A:=Matrix(n,n,t2,'shape'='symmetric');                       # table to Matrix
  LinearAlgebra:-Map(x->if x::list then x[2] else x end if,A); # keep only weights
  GraphTheory:-Graph(A);
end proc:

Download Graph.mw

I was working on this but see @Carl Love was first. Here's a slight variation.

Download partitions.mw

Some of the lambda's in vvalue were not formed with the double underline (use lprint vvalue to see this). After fixing this, simplify(...,size) executes in a reasonable time.
 

Download sol1det.mw

Hard to tell without an uploaded worksheet (use green up-arrow), but it looks like NN and MM weren't given values.

To evaluate y at a particular x value use eval (for evaluate). For example,

y:=0.1*x+0.3*x^2;
eval(y,x=2);

gives 1.4

To do a lot of them, set up y as a function (technically a procedure), and put the x values you want in a list:

y:=x->0.1*x+0.3*x^2;
xvals:=[seq(i,i=0..2,0.5)];
y~(xvals)

gives [0., 0.125, 0.400, 0.825, 1.400]

plottools:-exportplot can export .svg (should also do .tiff but I got an error message). The .svg file looks correct in CorelDraw.

Download exportplot.mw

From the ?dsolve,numeric help page "All IVP methods can be used for complex-valued IVPs with a real-valued independent variable". Not sure what it would mean for a complex independent variable - some sort of path through the complex plane would be needed.