restart;
d1:=3: d2:=3: d3:=5:
plots:-implicitplot3d(max(x^2+y^2+z^2-d1^2, (x-3)^2+y^2+z^2-d2^2, x^2+y^2+(z-4)^2-d3^2), x=0..4, y=-3..3, z=-2..3, style=surface, color=khaki, numpoints=500000, axes=normal, scaling=constrained);
               

The user lighting was used.

Addition. We can build everything together in one plot:

restart;
d1:=3: d2:=3: d3:=5:
A:=plots:-implicitplot3d(max(x^2+y^2+z^2-d1^2,(x-3)^2+y^2+z^2-d2^2,x^2+y^2+(z-4)^2-d3^2), x=0..4,y=-3..3,z=-2..3, style=surface, color=khaki, numpoints=500000, axes=normal, scaling=constrained):    
B:=plots:-implicitplot3d([x^2+y^2+z^2=d1^2,(x-3)^2+y^2+z^2=d2^2,x^2+y^2+(z-4)^2=d3^2],x=-5..6,y=-5..5,z=-4..9,color=["LightPink","LightGreen","LightBlue"], style= wireframe, thickness=0,  numpoints=5000):
plots:-display(A, B);
                

The  roots  returns roots over the field implied by the coefficients present (see help for details).  For example, if all the coefficients are rational, then the rational roots are computed. If you want all the roots (over the field of complex numbers) use  solve  command instead of  roots :

solve(x^3-x^2-8*x+8);

Matrix(3, (i,j)->M[i, j]+b[j] );

With this method, we explicitly indicate Maple the rule for obtaining elements of the new matrix.

Another method:
M := <1,2,3; 4,5,6; 7,8,9>;
b := <10 | 11 | 12>;
M1 := <b, b, b>;
M+M1;

Replace in your code  c__1^.5  by  sqrt(c__1)  or  c__1^(1/2) . The same for  c__2^.5

 :=  instead of  =  in the first line.

restart;
g:=x->(2*x-3)^``(-3);
g(x);
g:=x->(2*x-3)^`-3`;
g(x);
 

restart;
x := (-b + sqrt(b^2 - 4*a*c) ) /(2*a):
x := subsindets(x, sqrt, t->Z);
 

In Maple, you can do animation with only one parameter. Therefore,  r  should also depend on  theta. We obtain an animation of a space curve.

An example ( a  is the parameter of the animation):

r:=sqrt(theta):
plots:-animate(plots:-spacecurve, [[r*cos(theta),r*sin(theta),r*cos(theta)], theta=0..a, color=red, thickness=2], a=0..2*Pi);

In this example, we get a spiral lying in the plane  z=x .

Edit.

L:=`~`[int](convert(convert(series(x^x, x), polynom), list), x = 0 .. 1);
[L[1], seq(sign(L[i])*1/``(simplify(denom(abs(L[i]))^(1/i)))^i, i=2..nops(L))];
            

Of course, since we already know the pattern, we can make it simpler without referring to L:
 [1, seq((-1)^(n-1)*1/``(n)^n, n=2..6)];

Simply Maple is weak in solving some types of equations and in some plottings. Consider the following more interesting and less trivial example. Let  A  and  B  be two fixed points on a plane and construct the function  f  that associates with each point of the plane  P(x,y)  the sum of the distances from this point to points  A  and  B . Let  d  be the distance between points  A  and  B . Then the solution of the equation  f(x,y) = d  is the whole segment  [A, B] . Maple finds only the straight line passing through points A and B. With the plotting of this segment Maple also fails (due to the reason specified by Joe - no changes of signs):

A:=[1,1]:  B:=[5,4]:  d:=5: 
f:=(x,y)->sqrt((x-A[1])^2+(y-A[2])^2)+sqrt((x-B[1])^2+(y-B[2])^2):
solve(f(x,y)=5, {x,y});  # Incorrect solution. Should be 1<=y<=4
plots:-implicitplot(f(x,y)=d, x=0..6, y=0..6, gridrefine=5);  # No any curve
               

I explain to my students the differentiation of the function  x->a(x)^b(x)  using the identity
a(x)^b(x)=e^(b(x)*ln(a(x))  and then the rule of differentiation of the composition of functions.  Maple uses the same approch:

restart;
Student:-Calculus1:-ShowSolution(Diff(a(x)^b(x), x));

Since the result is expressed as an inert sum, it is best to use  value  command:

restart;
a:=unapply(rsolve({U(n)=U(n-1)+3*U(n-3), U(0) = 6, U(1) = 6, U(2) = 6, U(3)=24}, U(n)), n);
seq(value(a(n)), n=0..19);  # The first 20 terms
              

Download Digitalize_new.mw

The following code works:

restart;
ExY := map(convert,[WPKCPYWFYWCXHY, WPKCTEYFHCFEEE, HYCYTFHYYWWKCE, HYHHTFHKFEYHEH, HKXWYHFECTFFTF, HKWEHPPECWWHTC, HKXWYTPXHFHWYP, HKFEXCTFECXFKP, HKFEYPPEEEPHYW, HYTHCCFEWPPEXX, HKXKXCKFCHTEWK, HYEKPPCKFTWXXW], string);
Blist:=map(convert,[W,P,K,C,Y,F,X,H,T,E], string); # basis letters
nWords:=numelems(ExY);
Wletters:=length(ExY[1]);
Bletters:=numelems(Blist);
A:=Matrix(Bletters,Bletters):
with(ListTools):
for iW from 1 to nWords do 
  for iL from 1 to Wletters-1 do 
     fromLett:=ExY[iW][iL]; 
     toLett:=ExY[iW][iL+1]; 
     BfromLett:=Search(fromLett,Blist);
     BtoLett:=Search(toLett,Blist);
    A[BfromLett,BtoLett]:=A[BfromLett,BtoLett]+ 1; 
  end do; 
end do; 
A;

1. You have chosen an extremely inefficient way to create a list. Compare two ways:

 # The first way
t:=time():
L:=[$ 1..100000]:
time()-t;

# The second way
t:=time():  
 L := []: for i from 1 to 100000 do L := [op(L), i] end do:
time()-t;

2. I do not know how to refer to the labels programmatically. But if you assign any names to your equations, for example eq1, eq2, ... , eq10 , then 

[seq(eq||i, i=1..10)];