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

1. Conversion to  Matrix  can be easier (see the code below).

2. I agree that this piece "A particular solution of the system can be then written in the form   F(t) = S(t) C0 + P(t)  where  C0 is  n by 1 and  F(0) = C0 + P(0)." of Help is incorrect.

3. The constants  C1  and  C2  for the particular solution with the specified initial conditions  x0  can be obtained as follows (first we find the general solution  x ):

restart;
A:=Matrix(2, 2, [1, 4, 5, 1]);
x0:=Vector([2,1]); #initial condition
sol := DETools:-matrixDE(A, t);
M:=convert(sol[1], Matrix);  #convert to Matrix
x:=M.<C1,C2>;  # C1 and C2 are arbitrary constants
solve(eval(convert(x=~x0,list),t=0), [C1,C2])[];
simplify(eval(x, [%[],t=0])); # The check

dice := proc()
    local r;
    r:=rand(1..6);
    'r'()$2;
end proc:

Example of use:
seq([dice()], i=1..10);
             [6, 5], [3, 1], [5, 2], [3, 2], [2, 4], [3, 3], [1, 2], [5, 4], [5, 6], [2, 3]