I also did not understand exactly what you want to plot. You have a system of 5 unknown functions of the variable  t . Below I showed how you can plot the graph  V  vs  t. For other functions, everything is the same.

DEplot({diff(C(t), t) = rho*lambda*S(t)+rho*epsilon*lambda*V(t)+(1-q)*eta*I(t)-(mu+beta+chi)*C(t), diff(I(t), t) = (1-rho)*lambda*S(t)+(1-rho)*epsilon*lambda*V(t)+chi*C(t)-(mu+alpha+eta)*I(t), diff(R(t), t) = beta*C(t)+q*eta*I(t)-(mu+delta)*R(t), diff(S(t), t) = (1-p)*pi+phi*V(t)+delta*R(t)-(mu+lambda+vartheta)*S(t), diff(V(t), t) = p*pi+vartheta*S(t)-(epsilon*lambda+mu+phi)*V(t)}, {C(t), I(t), R(t), S(t), V(t)}, t = 0 .. 300, number = 5, stepsize = .1, [[S(0) = 8200, V(0) = 2800, C(0) = 200, I(0) = 210, R(0) = 200]], linecolour = t, axes = BOXED, scene = [t, V(t)]);

In fact, you study the monotonicity of the function  x-> (3*x-1)/(x^3+2*x+1)  in some intervals (for example for  x>0). For this instead  is  and  assume  command, I would just have used a plot, and to refine the intervals of increase and decrease, I would use the derivative:

restart;
f := x-> (3*x-1)/(x^3+2*x+1):
plot(f, 0..10, size=[700,300]);
R1:=solve({diff(f(x), x)>0, x>0});
R2:=solve({diff(f(x), x)<0, x>0});
evalf([R1, R2]);
 

I successfully restarted your worksheet and saved it again. Probably the cause of the problem was in the absence of any delimiters between the individual lines of your code (I had put them).

Plot_Diff_Temperature_new.mw

local D:
with(LinearAlgebra):
A:=RandomMatrix(1,1);
B:=RandomMatrix(1,5);
C:=RandomMatrix(5,1);
D:=RandomMatrix(5,5);
Matrix(2,2,{(1,1)=A,(1,2)=B,(2,1)=C,(2,2)=D});

Another way:

local D:
with(LinearAlgebra):
A:=RandomMatrix(1,1);
B:=RandomMatrix(1,5);
C:=RandomMatrix(5,1);
D:=RandomMatrix(5,5);
Matrix(2,2, ['A', 'B', 'C', 'D']):
M:=%;

Edit.

1. cutin  and  cutout  commands work only in 3D. If you try to run them in 2D, you receive an error message.

2. Yes,  project  command allows you to get shadows from some objects in 2D and 3D, but only with orthogonal projection. For an arbitrary projection, you can use  plottools[transform]  command.

with(plots): with(plottools):
P:=polygon([[0,0],[0,1],[1,1]], color=blue):
P1:=polygon([[0,0,1],[0,1,1],[1,1,1]], color=blue):
Q1:=polygon([[0,0,1],[0,1,1],[1,1,1]], color=grey):
display(cutin(P,1/2), axes=normal);  # An error
display(cutin(P1,1/2), axes=normal, view=[0..1,0..1,0..2]);
display(cutout(P,1/2), axes=normal);  # An error
display(cutout(P1,1/2), axes=normal, view=[0..1,0..1,0..2]);
display(cutout(P1,1/2), project(cutout(Q1,1/2), [[0,0,0],[1,0,0],[0,1,0]]), view=[0..1,0..1,0..2], axes=normal);  # Getting an shadow 

Formal parameters of the procedure: x  is a certain number, and  n  is a positive integer.

restart;
EVALF:=proc(x,n::posint)
local m;
m:=ilog10(evalf(x));
Digits:=m+n+3;
Float(round(Float(evalf(x), n)), -n);
end proc:
 

 
Examples of use:

EVALF(Pi, 3);
                                            3.142

EVALF(Pi, 15);
                                  3.141592653589793

Edit.

To obtain a numerical solution, you must specify the parameter values and add the initial conditions.

An example of the solution:

restart;
Sys:={-(9/14)*R(t)^(19/14)-(285/196)*(diff(R(t), t, t))/R(t)^(9/14)+(2565/2744)*(diff(R(t), t))^2/R(t)^(23/14) = -k*(4*lambda+1/a(t)^3),
R(t)= 6*((diff(a(t), t))^2/a(t)^2+(diff(a(t),t,t))/a(t))}:
k:=2: lambda:=1:
Sol:=dsolve(Sys union {R(0)=0.5,a(0)=1.5,D(R)(0)=0,D(a)(0)=1}, {R(t),a(t)}, numeric);
plots:-odeplot(Sol,[[t,R(t)],[t,a(t)]], t=0..7, color=[red,blue]);

Another way to solve the problem is to call  RealDomain  package. Then all calculations Maple does in the real domain:

with(RealDomain):
P:=plot([(x^3-x^2)^(1/3), x-1/3, [[0,0]], [[1,0]], [[2/3,0]], [[0,-4^(1/3)/3]], [t,-4^(1/3)/3,t=0..2/3], [2/3,t,t=0..-4^(1/3)/3] ], x=-2.5..2.5, y=-2.3..2.3, style=[line$2,point$4,line$2], linestyle=[1,3,2$2], color=[red,black,red$4,black$2], thickness=[2,1], symbol=solidcircle, symbolsize=[8,8,11,11], size=[600,500], scaling=constrained, caption = typeset("A plot of   ", ``(x^3-x^2)^(1/3)), captionfont=[times,16]):
T:=plots:-textplot([[0.75,0.15,2/3], [0.15,-0.7,-4^(1/3)/3], [0.9,1,y=x-1/3]]):
plots:-display(P,T);

This can be done in several ways. Here are two ways:

convert(3^665, string)[-10..-1];  # The last 10 digits
                                                               "7817543843"

seq(convert(3^665, base, 10)[i], i=10..1,-1);
                                                         7, 8, 1, 7, 5, 4, 3, 8, 4, 3

Try this:

Vv := Eigenvalues(prod);
simplify(Vv);

Input:
 

SymmetricReduction[
 2*i1^4*i2*i3 + 2*i1^4*i2*i4 + 2*i1^4*i2*i5 + 2*i1^4*i3*i4 + 
  2*i1^4*i3*i5 + 2*i1^4*i4*i5 + 4*i1^3*i2^2*i3 + 4*i1^3*i2^2*i4 + 
  4*i1^3*i2^2*i5 + 4*i1^3*i2*i3^2 + 6*i1^3*i2*i3*i4 + 
  6*i1^3*i2*i3*i5 + 4*i1^3*i2*i4^2 + 6*i1^3*i2*i4*i5 + 
  4*i1^3*i2*i5^2 + 4*i1^3*i3^2*i4 + 4*i1^3*i3^2*i5 + 4*i1^3*i3*i4^2 + 
  6*i1^3*i3*i4*i5 + 4*i1^3*i3*i5^2 + 4*i1^3*i4^2*i5 + 
  4*i1^3*i4*i5^2 + 2*i1^2*i2^3*i3 + 2*i1^2*i2^3*i4 + 2*i1^2*i2^3*i5 + 
  4*i1^2*i2^2*i3^2 + 6*i1^2*i2^2*i3*i4 + 6*i1^2*i2^2*i3*i5 + 
  4*i1^2*i2^2*i4^2 + 6*i1^2*i2^2*i4*i5 + 4*i1^2*i2^2*i5^2 + 
  2*i1^2*i2*i3^3 + 6*i1^2*i2*i3^2*i4 + 6*i1^2*i2*i3^2*i5 + 
  6*i1^2*i2*i3*i4^2 + 24*i1^2*i2*i3*i4*i5 + 6*i1^2*i2*i3*i5^2 + 
  2*i1^2*i2*i4^3 + 6*i1^2*i2*i4^2*i5 + 6*i1^2*i2*i4*i5^2 + 
  2*i1^2*i2*i5^3 + 2*i1^2*i3^3*i4 + 2*i1^2*i3^3*i5 + 
  4*i1^2*i3^2*i4^2 + 6*i1^2*i3^2*i4*i5 + 4*i1^2*i3^2*i5^2 + 
  2*i1^2*i3*i4^3 + 6*i1^2*i3*i4^2*i5 + 6*i1^2*i3*i4*i5^2 + 
  2*i1^2*i3*i5^3 + 2*i1^2*i4^3*i5 + 4*i1^2*i4^2*i5^2 + 
  2*i1^2*i4*i5^3 + 2*i1*i2^3*i3*i4 + 2*i1*i2^3*i3*i5 + 
  2*i1*i2^3*i4*i5 + 4*i1*i2^2*i3^2*i4 + 4*i1*i2^2*i3^2*i5 + 
  4*i1*i2^2*i3*i4^2 + 6*i1*i2^2*i3*i4*i5 + 4*i1*i2^2*i3*i5^2 + 
  4*i1*i2^2*i4^2*i5 + 4*i1*i2^2*i4*i5^2 + 2*i1*i2*i3^3*i4 + 
  2*i1*i2*i3^3*i5 + 4*i1*i2*i3^2*i4^2 + 6*i1*i2*i3^2*i4*i5 + 
  4*i1*i2*i3^2*i5^2 + 2*i1*i2*i3*i4^3 + 6*i1*i2*i3*i4^2*i5 + 
  6*i1*i2*i3*i4*i5^2 + 2*i1*i2*i3*i5^3 + 2*i1*i2*i4^3*i5 + 
  4*i1*i2*i4^2*i5^2 + 2*i1*i2*i4*i5^3 + 2*i1*i3^3*i4*i5 + 
  4*i1*i3^2*i4^2*i5 + 4*i1*i3^2*i4*i5^2 + 2*i1*i3*i4^3*i5 + 
  4*i1*i3*i4^2*i5^2 + 2*i1*i3*i4*i5^3 + 4*i1^3*i2*i3 + 4*i1^3*i2*i4 + 
  4*i1^3*i2*i5 + 4*i1^3*i3*i4 + 4*i1^3*i3*i5 + 4*i1^3*i4*i5 + 
  8*i1^2*i2^2*i3 + 8*i1^2*i2^2*i4 + 8*i1^2*i2^2*i5 + 8*i1^2*i2*i3^2 + 
  12*i1^2*i2*i3*i4 + 12*i1^2*i2*i3*i5 + 8*i1^2*i2*i4^2 + 
  12*i1^2*i2*i4*i5 + 8*i1^2*i2*i5^2 + 8*i1^2*i3^2*i4 + 
  8*i1^2*i3^2*i5 + 8*i1^2*i3*i4^2 + 12*i1^2*i3*i4*i5 + 
  8*i1^2*i3*i5^2 + 8*i1^2*i4^2*i5 + 8*i1^2*i4*i5^2 + 4*i1*i2^3*i3 + 
  4*i1*i2^3*i4 + 4*i1*i2^3*i5 + 8*i1*i2^2*i3^2 + 12*i1*i2^2*i3*i4 + 
  12*i1*i2^2*i3*i5 + 8*i1*i2^2*i4^2 + 12*i1*i2^2*i4*i5 + 
  8*i1*i2^2*i5^2 + 4*i1*i2*i3^3 + 12*i1*i2*i3^2*i4 + 
  12*i1*i2*i3^2*i5 + 12*i1*i2*i3*i4^2 + 48*i1*i2*i3*i4*i5 + 
  12*i1*i2*i3*i5^2 + 4*i1*i2*i4^3 + 12*i1*i2*i4^2*i5 + 
  12*i1*i2*i4*i5^2 + 4*i1*i2*i5^3 + 4*i1*i3^3*i4 + 4*i1*i3^3*i5 + 
  8*i1*i3^2*i4^2 + 12*i1*i3^2*i4*i5 + 8*i1*i3^2*i5^2 + 4*i1*i3*i4^3 + 
  12*i1*i3*i4^2*i5 + 12*i1*i3*i4*i5^2 + 4*i1*i3*i5^3 + 4*i1*i4^3*i5 + 
  8*i1*i4^2*i5^2 + 4*i1*i4*i5^3 + 4*i2^3*i3*i4 + 4*i2^3*i3*i5 + 
  4*i2^3*i4*i5 + 8*i2^2*i3^2*i4 + 8*i2^2*i3^2*i5 + 8*i2^2*i3*i4^2 + 
  12*i2^2*i3*i4*i5 + 8*i2^2*i3*i5^2 + 8*i2^2*i4^2*i5 + 
  8*i2^2*i4*i5^2 + 4*i2*i3^3*i4 + 4*i2*i3^3*i5 + 8*i2*i3^2*i4^2 + 
  12*i2*i3^2*i4*i5 + 8*i2*i3^2*i5^2 + 4*i2*i3*i4^3 + 
  12*i2*i3*i4^2*i5 + 12*i2*i3*i4*i5^2 + 4*i2*i3*i5^3 + 4*i2*i4^3*i5 + 
  8*i2*i4^2*i5^2 + 4*i2*i4*i5^3 + 4*i3^3*i4*i5 + 8*i3^2*i4^2*i5 + 
  8*i3^2*i4*i5^2 + 4*i3*i4^3*i5 + 8*i3*i4^2*i5^2 + 4*i3*i4*i5^3 + 
  i1^4 + 3*i1^3*i2 + 3*i1^3*i3 + 3*i1^3*i4 + 3*i1^3*i5 + 
  3*i1^2*i2^2 + 6*i1^2*i2*i3 + 6*i1^2*i2*i4 + 6*i1^2*i2*i5 + 
  3*i1^2*i3^2 + 6*i1^2*i3*i4 + 6*i1^2*i3*i5 + 3*i1^2*i4^2 + 
  6*i1^2*i4*i5 + 3*i1^2*i5^2 + i1*i2^3 + 3*i1*i2^2*i3 + 
  3*i1*i2^2*i4 + 3*i1*i2^2*i5 + 3*i1*i2*i3^2 + 10*i1*i2*i3*i4 + 
  10*i1*i2*i3*i5 + 3*i1*i2*i4^2 + 10*i1*i2*i4*i5 + 3*i1*i2*i5^2 + 
  i1*i3^3 + 3*i1*i3^2*i4 + 3*i1*i3^2*i5 + 3*i1*i3*i4^2 + 
  10*i1*i3*i4*i5 + 3*i1*i3*i5^2 + i1*i4^3 + 3*i1*i4^2*i5 + 
  3*i1*i4*i5^2 + i1*i5^3 + 4*i2^2*i3*i4 + 4*i2^2*i3*i5 + 
  4*i2^2*i4*i5 + 4*i2*i3^2*i4 + 4*i2*i3^2*i5 + 4*i2*i3*i4^2 + 
  4*i2*i3*i5^2 + 4*i2*i4^2*i5 + 4*i2*i4*i5^2 + 4*i3^2*i4*i5 + 
  4*i3*i4^2*i5 + 4*i3*i4*i5^2, {i2, i3, i4, i5}]

Output:

{i1^4 + (48 i1 + 24 i1^2 + 12 (-4 i1 - 2 i1^2) + 
     6 (-8 i1 - 4 i1^2 - 2 (-4 i1 - 2 i1^2)) + 
     4 (-4 - 12 i1 - 6 i1^2 - 5 (-4 i1 - 2 i1^2) - 
        2 (-8 i1 - 4 i1^2 - 2 (-4 i1 - 2 i1^2)))) i2 i3 i4 i5 + 
  3 i1^3 (i2 + i3 + i4 + i5) + (12 + 3 (-8 - 2 (-4 - 2 i1) - 4 i1) + 
     7 (-4 - 2 i1) + 6 i1) i2 i3 i4 i5 (i2 + i3 + i4 + i5) + 
  3 i1^2 (i2 + i3 + i4 + i5)^2 + 
  i1 (i2 + i3 + i4 + i5)^3 + (4 i1^3 + 2 i1^4) (i2 i3 + i2 i4 + 
     i3 i4 + i2 i5 + i3 i5 + i4 i5) + (8 i1^2 + 4 i1^3) (i2 + i3 + 
     i4 + i5) (i2 i3 + i2 i4 + i3 i4 + i2 i5 + i3 i5 + 
     i4 i5) + (4 i1 + 2 i1^2) (i2 + i3 + i4 + i5)^2 (i2 i3 + i2 i4 + 
     i3 i4 + i2 i5 + i3 i5 + i4 i5) + (8 i1 + 4 i1^2 + 
     2 (-4 i1 - 2 i1^2)) (i2 i3 + i2 i4 + i3 i4 + i2 i5 + i3 i5 + 
     i4 i5)^2 + (4 i1 + 12 i1^2 + 6 i1^3 + 
     3 (-8 i1^2 - 4 i1^3)) (i2 i3 i4 + i2 i3 i5 + i2 i4 i5 + 
     i3 i4 i5) + (4 + 12 i1 + 6 i1^2 + 5 (-4 i1 - 2 i1^2) + 
     2 (-8 i1 - 4 i1^2 - 2 (-4 i1 - 2 i1^2))) (i2 + i3 + i4 + 
     i5) (i2 i3 i4 + i2 i3 i5 + i2 i4 i5 + i3 i4 i5) + (4 + 
     2 i1) (i2 + i3 + i4 + i5)^2 (i2 i3 i4 + i2 i3 i5 + i2 i4 i5 + 
     i3 i4 i5) + (8 + 2 (-4 - 2 i1) + 4 i1) (i2 i3 + i2 i4 + i3 i4 + 
     i2 i5 + i3 i5 + i4 i5) (i2 i3 i4 + i2 i3 i5 + i2 i4 i5 + 
     i3 i4 i5), 0}

Check in Maple:

P:=i1^4 + (48*i1 + 24*i1^2 + 12*(-4*i1 - 2*i1^2) + 
     6*(-8*i1 - 4*i1^2 - 2*(-4*i1 - 2*i1^2)) + 
     4*(-4 - 12*i1 - 6*i1^2 - 5*(-4*i1 - 2*i1^2) - 
        2*(-8*i1 - 4*i1^2 - 2*(-4*i1 - 2*i1^2))))*i2*i3*i4*i5 + 
  3*i1^3*(i2 + i3 + i4 + i5) + (12 + 3*(-8 - 2*(-4 - 2*i1) - 4*i1) + 
     7*(-4 - 2*i1) + 6*i1)*i2*i3*i4*i5*(i2 + i3 + i4 + i5) + 
  3*i1^2*(i2 + i3 + i4 + i5)^2 + 
  i1*(i2 + i3 + i4 + i5)^3 + (4*i1^3 + 2*i1^4)*(i2*i3 + i2*i4 + 
     i3*i4 + i2*i5 + i3*i5 + i4*i5) + (8*i1^2 + 4*i1^3)*(i2 + i3 + 
     i4 + i5)*(i2*i3 + i2*i4 + i3*i4 + i2*i5 + i3*i5 + 
     i4*i5) + (4*i1 + 2*i1^2)*(i2 + i3 + i4 + i5)^2*(i2*i3 + i2*i4 + 
     i3*i4 + i2*i5 + i3*i5 + i4*i5) + (8*i1 + 4*i1^2 + 
     2*(-4*i1 - 2*i1^2))*(i2*i3 + i2*i4 + i3*i4 + i2*i5 + i3*i5 + 
     i4*i5)^2 + (4*i1 + 12*i1^2 + 6*i1^3 + 
     3*(-8*i1^2 - 4*i1^3))*(i2*i3*i4 + i2*i3*i5 + i2*i4*i5 + 
     i3*i4*i5) + (4 + 12*i1 + 6*i1^2 + 5*(-4*i1 - 2*i1^2) + 
     2*(-8*i1 - 4*i1^2 - 2*(-4*i1 - 2*i1^2)))*(i2 + i3 + i4 + 
     i5)*(i2*i3*i4 + i2*i3*i5 + i2*i4*i5 + i3*i4*i5) + (4 + 
     2*i1)*(i2 + i3 + i4 + i5)^2*(i2*i3*i4 + i2*i3*i5 + i2*i4*i5 + 
     i3*i4*i5) + (8 + 2*(-4 - 2*i1) + 4*i1)*(i2*i3 + i2*i4 + i3*i4 + 
     i2*i5 + i3*i5 + i4*i5)*(i2*i3*i4 + i2*i3*i5 + i2*i4*i5 + 
     i3*i4*i5);

expand(P) - Teller;
length(Teller);
length(P);

                                              0
                                           4308
                                            969
 

All the plotting for the first curve. For the second curve everything is done similarly.

f:=x->sqrt(abs(x-1));
y:=f(a)+D(f)(a)*(x-a);
plots:-animate(plot,[ [y, [[a,f(a)]] ], x=a-2..a+2, -1..3, style=[line,point], color=red, thickness=2, symbolsize=15, scaling=constrained], a=3..1.00001, background=plot(f, -1..5, color=blue), frames=90);

The  Sol  procedure for Newton method.
Formal parameters: eta  (should be >1) and eps (positive number specifying an accuracy).

restart;
Sol:=proc(eta,eps::positive)
local y, n, k;
global F;
k:=-ilog10(eps);
Digits:=k+2;
if eta<=1 then error "Should be eta>1" fi;
F:=(y,eta)->eta*(2+9*exp(-eta*y))-22*eta*y;
y[0]:=0.;
for n from 1 do
y[n]:=y[n-1]-F(y[n-1],eta)/D[1](F)(y[n-1],eta);
if abs(y[n]-y[n-1])<eps then break fi;
od;
evalf[max(k,1)](y[n]);
end proc:

Example of use and check:

y0:=Sol(2, 0.0001);
y01:=fsolve(F(y,2), y=0..1);
abs(y0-y01);
plot(F(y,2), y=0..1);

The  Bis  procedure for the bisection method.

restart;
Bis:=proc(eta,eps::positive)
local a, b, c, n, k;
global F;
k:=-ilog10(eps);
Digits:=k+2;
if eta<=1 then error "Should be eta>1" fi;
F:=(y,eta)->eta*(2+9*exp(-eta*y))-22*eta*y;
a:=0.;  b:=1.;
for n from 1 do
c:=(a+b)/2;
if F(a,eta)*F(c,eta)<0 then b:=c else a:=c fi;
if abs(a-b)<eps then break fi;
od;
evalf[max(k,1)](c);
end proc:

Example of use:

y0:=Bis(2, 0.0001);

Edit.

You do not need  isolve  command. (j-1)/m<=t<j/m  is equivalent to  (j-1)<=t*m<j , therefore
j-1=floor(t*m) ,   j=floor(t*m) + 1

Addition. If you still want to use  isolve  command, then it's easy to do as some procedure:

F:=(m,t)->isolve({t>=(j-1)/m, t<j/m}):

Example of use:

F(2.5, 0.7);
                                          {j = 2}

You have a second-order equation, so for each solution you need two initial conditions:

ode1 := diff(x(t), t, t) = (5*9.80665)*sin((1/6)*Pi)-(10*(10-sqrt(x(t)^2+25)))*x(t)/sqrt(x(t)^2+25)-diff(x(t),t);
DEtools:-DEplot(ode1, x(t), t = -2 .. 2, x = -8 .. 8, color = blue, stepsize = 0.5e-1, [seq([x(0) = (1/4)*k, D(x)(0)=0], k = -20 .. 20, 4)], linecolor=red);

Arrows are plotted only for first-order equations.