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)];

From the link   webpage  you specified, you missed the beginning of the code, in which all parameters  (a, b, h, n, f)  were set. Here is the full code for your version (you have changed only the parameter  n ):

Download Adjust.mw

The equation  DEF_Profit1  in addition to the unknown  S  contains 3 parameters  beta, s, C. If you specify the values of these parameters, you can easily find all the roots, for example:

sol1 := solve(eval(DEF_Profit1, [beta = 3, s = 0.1, C = 2]), S);
                        sol1 := 0.004990229958, 29.61213284, -3.574031930, -26.04309114

By default, axes=normal  option works. So just remove  axes=box  option from your code:

plot(sin(x), x=0..2*Pi);

Addition.
 If you want the boxed effect and that the tickmarks are on the line  y=0 , then do so:

plots:-display(plot(sin(x), x=0..2*Pi), plot([[0,-1],[2*Pi,-1],[2*Pi,1],[0,1]], color=black));
                

Edit.         

If you do not have additional information about the parameters  chi, lambda1, lambda2, phi, psi, rho, you can simply use  simplify  command:

with(LinearAlgebra):

A:=Matrix([[phi,conjugate(psi),chi,conjugate(rho)],[psi,-conjugate(phi),rho,-conjugate(chi)],[lambda1*phi,conjugate(lambda1)*conjugate(psi),lambda2*chi,conjugate(lambda2)*conjugate(rho)],[lambda1*psi,-conjugate(lambda1)*conjugate(phi),lambda2*rho,-conjugate(lambda2)*conjugate(chi)]]);
indets(A, name);
Determinant(A);
simplify(%);
                           

Try  colorscheme  option for curves. An example:

restart;
plots:-spacecurve([cos(t), sin(t), 0.3*t], t=0..4*Pi, colorscheme=[blue, green], thickness=3, axes=normal);
                            

Unfortunately, this does not work when solving differential equations with numeric option (plotting by  plots:-odeplot  command). In this case, use a workaround as below (the same curve as above):

restart;
sol:=dsolve({diff(x(t),t)=-sin(t), diff(y(t),t)=cos(t),diff(z(t),t)=0.3, x(0)=1,y(0)=0,z(0)=0}, numeric):
plots:-pointplot3d([seq(eval([x(t),y(t),z(t)],sol(i)), i = 0..12.56,0.2)], colorscheme = [blue, green], style = line, thickness=3, axes=normal);
 

You can specify  g(x)  as a piecewise-function and then use in any calculations, for example:

f(x):=x+1;
g(x):=piecewise(x<>2,(x+1)*(x-2)/(x-2), x=2,undefined);
simplify(f(x)*(g(x)+1));
limit(g(x), x=2);
diff(g(x), x);
int(g(x), x=1..3);
discont(g(x), x);
                                   

With  for  loop:

restart;
r:=rand(1..100):
k:=0:
for i from 1 to 20 do
L[i]:=r();
if type(L[i], even) then k:=k+1; L1[k]:=L[i] fi;
od:
L:=convert(L, list);
L1:=convert(L1, list);
`+`(L1[])/nops(L1);

It is easier to do this in another way:

restart;
r:=rand(1..100):
L:=[seq(r(), i=1..20)];
L1:=select(t->type(t, even), L);
`+`(L1[])/nops(L1);

If you want new random numbers to be generated with each call, then insert the  line  randomize():  after restart;

Programmatically, I found only a method using  plots:-textplot  command:

plots:-textplot([[-0.1,0, max,font=[times,bold,20]], [0.1,0, f(x),font=[times,20]], [0,-0.1, x in [2,3]]], size=[140,60], axes=none);