Here is the procedure that inserts a rational number with the lowest denominator (vv's assumption) between any 2 list items. If there are several such numbers, then the smallest number is inserted.

FindSamples:=proc(sourcesamples)
local N, P;
N:=nops(sourcesamples);
P:=proc(a,b)
local a1, b1, m1, n, m;
if a=b then error "Should be a<>b" fi;
a1,b1:=op(convert(sort([a,b],(x,y)->evalf(x)<evalf(y)),rational));
for n from 1 do
m1:=a1*n;
m:=`if`(type(m1,integer),m1+1,ceil(m1));
if is(m/n>a1) and is(m/n<b1) then return m/n fi;
od;
end proc:
[ceil(sourcesamples[1])-1, seq(op([sourcesamples[i],P(sourcesamples[i],sourcesamples[i+1])]), i=1..N-1),sourcesamples[N],floor(sourcesamples[N])+1];
end proc:

The usage for your example:
sourcesamples :=[evalf(-1-sqrt(7)), -2, 1, evalf(-1 + sqrt(7)), 2];
FindSamples(sourcesamples);
                    

The main role in this procedure is played by the subprocedure  P , which returns a rational number between  a  and  b  with the lowest denominator. Here is a specific example of the application of this procedure: Take 2 numbers 3.14 and 3.15 (these are approximations of the number  Pi  with 3 digits with a deficiency and an excess) and find the intermediate number with the lowest denominator:

P:=proc(a,b) 
local a1, b1, m1, n, m; 
if a=b then error "Should be a<>b" fi;
a1,b1:=op(convert(sort([a,b],(x,y)->evalf(x)<evalf(y)),rational)); 
for n from 1 do 
m1:=a1*n; 
m:=`if`(type(m1,integer),m1+1,ceil(m1)); 
if is(m/n>a1) and is(m/n<b1) then return m/n fi; 
od;
end proc:
P(3.14, 3.15);

                         
 
Edit. The procedure  P  has been improved. Now it works much faster.

First, we complete the squares to estimate the ranges for the loop variables.
There are 2 integer solutions:

restart;
Eq:=7*(x^2+y^2+z^2)+2*x+4*y-8*z-8;
Student:-Precalculus:-CompleteSquare(Eq, [x,y,z]);
Eq1:=%/7; k:=0: T:=table(): 
for x from -2 to 2 do
for y from -2 to 2 do
for z from -2 to 2 do
if Eq=0 then k:=k+1; T[k]:=[x,y,z] fi;
od: od: od:
T:=convert(T,list);

Two syntax errors:
1. N/2  should be a positive integer.
2. Should be  add  instead of  sum.

Your example:

N:=10:
A:=Array(1..N/2);
add(A);

Here is your example automatically. It works for an arbitrary list:

restart;
L:= [-1+sqrt(7), -1 - sqrt(7)]:

L1:=sort(L,(x,y)->evalf(x)<=evalf(y));
m:=floor(L1[1]);  M:=ceil(L1[-1]);
[$ m..M, op(L)];
samples:=sort(%, (x,y)->evalf(x)<=evalf(y));
evalf(samples);

Edit. If you need to insert something between the list items, you can do so:

[seq(op([A,s]), s=samples), A];

solve(a*x^2 + b*x + c=0, x, parametric=true, real);

Here is a simple procedure that allows you to beautifully plot graphs in 3D. Required parameters: Names - the list of  vertex names, Coords - the list of vertex coordinates, Edges - the set of edges. Opt  is an optional parameter that determines the graphic properties of objects (the color of the vertices, the size of the vertices (as a fraction of the average distance between the vertices), the color of the edges, the thickness of the edges).

 

Download DrawGraph3D_2.mw

Animation of this picture:

plots:-animate(plots:-display,[plottools:-rotate(%,phi,[[0,0,0],[0,0,1]])], phi=0..2*Pi, frames=60);

Edit.

Use the  dimension=3  option.

Example:

with(GraphTheory):
H := CompleteGraph(4);
DrawGraph(H, style = spring, dimension = 3);

plots:-animate(plot3d, [[x, (cos(x)+1)*cos(phi), (cos(x)+1)*sin(phi)], x = 0 .. 4*Pi, phi = (1/2)*Pi .. t, style = surface, color = khaki, labels = [x, y, z], lightmodel = light1], t = (1/2)*Pi .. (-3*Pi)*(1/2), frames = 60, axes = normal, view = [-1 .. 13.7, -2.4 .. 2.4, -2.4 .. 2.4])

I guess that OP only meant to animate the process of painting one area (the curves themselves should not move):

restart; 
with(plots): 
animate(shadebetween, [sin(x), cos(x), x = (1/4)*Pi .. t, positiveonly], 
t = (1/4)*Pi .. 5*Pi*(1/4),background = plot([sin(x), cos(x)], x = 0 .. 2*Pi, 
gridlines = true, thickness = 3, labels = ["", ""]), scaling = constrained, size = [900, 300])

The answer is obvious - the function  2*sqrt(xyz)  does not depend on  x . xyz is just a symbol unrelated to x .
In Maple, this can be tested using the  depends  function:

depends(x^2,x);
depends(xyz,x);

                                     true
                                    false

In the usual sense, it is impossible to plot a function of three variables, since it requires 3 dimensions for the arguments and one dimension for the values of the function, that is, we need 4-dimensional space. I propose the following way to solve problem: we fix 1 variable, for example  z , and then plot the function of two variables  x  and  y . Using  Explore command we can trace how this plot changes as  z  changes.

The last 2 lines of code show that for some values of arguments your function takes imaginary values. For such values, the plot is not built.

Download Explore_plot3d.mw

To solve your problem, I used the simplest formulas for numerical differentiation. Your expression  diff(u(y, t), y)+diff(diff(u(y, t), y), t)  is calculated for the 4 indicated values  b   by replacing the derivatives at the point  (y,t) = (0,1)  with the corresponding difference ratios with the step  h=0.001 : 

Download ND.mw

Example:

combinat:-permute([a,b,c,d]);

       [[a, b, c, d], [a, b, d, c], [a, c, b, d], [a, c, d, b], [a, d, b, c], [a, d, c, b], [b, a, c, d], [b, a, d, c], [b, c, a, d], [b, c, d, a], [b, d, a, c], [b, d, c, a], [c, a, b, d], [c, a, d, b], [c, b, a, d], [c, b, d, a], [c, d, a, b], [c, d, b, a], [d, a, b, c], [d, a, c, b], [d, b, a, c], [d, b, c, a], [d, c, a, b], [d, c, b, a]] 

In fact, there are many more solutions than indicated by vv, for example

restart;
sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0];
fsolve(eval(sys,[a=10,b=10.4]), {x=0..infinity,y=0..infinity,z=0..infinity});
;

                      sys := [x^2-x*y+y^2-1 = 0, -a^2+y^2-y*z+z^2 = 0, -b^2+x^2-x*z+z^2 = 0]
                                    {x =0.1989730189, y = 1.084528287, z = 10.49805888}

Further, we strictly show that there are positive solutions  (x, y, z)  for arbitrary a=b>=1. From considerations of continuity, this implies that for any number  a>=1 and a sufficiently small positive number c  for  (a, b=a+c)  there are also a positive solutions  (x, y, z) .

restart;
sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0];
Sol:=simplify(solve(eval(sys,[b=a]), [x,y,z], explicit));
a in solve(a>=1 and `or`(seq((eval(x,k)>0 and eval(y,k)>0 and eval(z,k)>0), k=%)));

 The result:                                          a in RealRange(1, Open(infinity))

Numerical experiments show that all solutions seem to be in the strip  1<=a<=b<a+1/2 . This does not mean that every point of this strip is a solution. This is only a necessary condition.

restart;
a := 1: b:= 2: l:= 1: alpha[1]:= 1: alpha[2]:= 3:
syspde:= [diff(u(x, t), t)-a+u(x, t)-u(x, t)^2*v(x, t)-alpha[1]*(diff(u(x, t), x$2)) = 0, diff(v(x, t), t)-b+u(x, t)^2*v(x, t)-alpha[2]*(diff(v(x, t), x$2)) = 0]:
Bcs:= [u(x,0)=1,v(x,0)=1,D[1](u)(-l, t) = 0,D[1](u)(l, t) = 0,D[1](v)(-l, t) = 0,D[1](v)(l, t) = 0]:
sol:=pdsolve(syspde, Bcs, numeric);

sol:-value(u(x, t)+diff(u(x, t), x)+diff(u(x, t), t), t=0.7)(0.5);

                                          sol:=module() ... end module
               [x = 0.5, t = 0.7, u(x, t)+diff(u(x, t), x)+diff(u(x, t), t) = 2.34305224942915]