It seems that HeatMap cannot be used with display,insequence.
[Also, the PLOT structure generated by HeatMap  contains some undocumented stuff].

But Explore works; however sometimes it may happen that the gridlines are missing (I don't know why).

restart;
with(Statistics):# with(LinearAlgebra): randomize(): with(plots):
f:= rand( 0..1 ):
ee:=proc()
  local RM := Matrix(10, 10, (i, j) ->2*f()-1 );
  Threads:-Sleep(0.5);
  HeatMap(RM, axis = [gridlines = [10, color = blue], thickness = 2]) ;
end:
Explore(ee(n), parameters=[[n=1..10, animate]], numframes=10);

1. You have many parameters which must be given numeric values in order to obtain a numeric solution.

indets([sys],name(Not(constant))) minus {t};

2. Please check if you mean pi or Pi.

3. Use
sol := dsolve(...)
because p is already a parameter.

4. You also have a syntax error in odeplot.

plots[odeplot](sol, [[t, S(t)], [t, V(t)], [t, C(t)], [t, I(t)], [t, R(t)]], t = 0 .. 200, numpoints = 25)

(do you really want numpoints=25?).

For such problems you should use one of the packages combinat or Iterator.

P.S. O simple method to test e.g. that i,j,k,l,m,n,o,p are distinct is to use nops({i,j,k,l,m,n,o,p}) = 8.

pdsolve looks for classical solutions, C^1 in this case. Your pde has not such solutions: the solution without bc, ic is u(x,t)=_F(t-x).

The solution given by Mathematica is not differentiable.
 

f:=int(sqrt(sin(x)),x):
a:=limit(f,x=Pi/2,left):
answer:=piecewise(x<Pi/2,f-a, x=Pi/2,0, -eval(f-a,x=Pi-x));

n:=4:
Z:=Matrix(n):  # zero matrix
LinearAlgebra:-MatrixExponential(Z);

This should be mathematically obvious because exp(Z) = I + Z^1/1! + Z^2/2! + ... = I

You should be aware that an expression has many equivalent forms.

SUBS:=proc(f,x)
local F,N:=0, g;
g:=subs(x=F(),f);
F:=proc() N:=N+1; x[N] end;
eval(g)
end:

SUBS(exp(-I*Pi*(n+2*n*(m-1))/m), n);

Probably you want integer (or rational) points, otherwise it's trivial.
Such points do not exist. Use isolve.

assume has limited capabilities especially for more than one variables.

restart;
f:=n -> (3*n-1)/(n^3+2*n+1):
is(f(n)<f(m)) assuming n>1,m>n;                  #false,OK
is(f(n)<f(m)) assuming n>1,m>n,n::integer;       #false,OK
is(f(n)>f(m)) assuming n>1,m>n;                  #FAIL,should be false
is(f(n)>f(m)) assuming n>1,m>n,n::integer;       #FAIL,should be true
g:=D(f):
is(g(n) < 0) assuming n>1;                 #false,OK
is(g(n) < 0) assuming n>1, n::integer;     #true,OK
is(g(n) < 0) assuming n>=2;                #true,OK

The value of each double integral is actually +oo.

The numerical routines may have problems in such cases. [Here one of the factors is small and the other is large near 0].

A0:=Matrix(1,1,(i,j)->a):
B0:=Matrix(1,5,(i,j)->b):
C0:=Matrix(5,1,(i,j)->c):
D0:=Matrix(5,5,(i,j)->d):
<<A0|B0>,<C0|D0>>;

With your definitions of the exports, :-Proc1   etc  must be defined after the module.
Actually this is equivalent to defining  export Proc1 := :-Proc1    etc

So, this works:

#ff:=x->2*x+3;
M := module()
  option package;
  export f;
  f:=eval(:-ff);
end module:
ff:=x->2*x+3:
with(M):
f(x): # 2*x+3
unwith(M):
#f(x);
ff:=x -> 5*x:
#f(x);
with(M):
f(x);  #5*x

But if you uncomment the fist line, without a restart f is not changed.

(Corrected)

3 &^ 665 mod 10^10;
       7817543843

This also works for huge numbers:

3 &^ 123456789 mod 10^10;
       5225122483