I assume that you want to do this so that you can later recall the alpha's and beta's.  That is, the older versions of alpha and beta don't appear to be necessary in the algorithm.  If that is the case, then you could store them in a table, using a global counter to increment the index.  Say

LocSubAlg := module()
global cnt, Alpha, Beta; 
...
   someproc := proc()
       alpha := ... ;
       beta := ... ;
       cnt := cnt+1;
       Alpha[cnt] := alpha;
       Beta[cnt] := beta;
    end proc:
end module:  

Actually I'd probably do that a bit differently, making Alpha, Beta and cnt module local and then returning them as the final step (in Start). But this should get you started.

Here's an efficiency suggestion.  While I don't know what type of input you expect to feed this, on the few tests I tried, the algorithm generated a lot of Matrices.  You can make this much more efficient by precomputing and reusing Matrices.  For example, consider the assignment to alpha1 in step2:

        alpha1 := ( Matrix([[1, 0, 0],
                            [0, 0, 1],
                            [0, -1, 0]])
                    . Matrix([[1, 0, 0],
                              [0, 1, ((-1)*a1*(B[2,2]))],
                              [0, 0, 1]])
                    . Matrix([[1, 0, 0],
                              [(B[2,1])*(B[2,2]), 1, 0],
                              [0, 0, 1]])
                  );

You don't need to compute that Matrix product with each iteration.  Instead, compute that symbolically to get

                         [       1             0        0     ]
                         [                                    ]
               alpha1 := [       0             0        1     ]
                         [                                    ]
                         [-B[2, 1] B[2, 2]    -1    a1 B[2, 2]]

So you could do

     alpha1 := Matrix(3,3,[1,0,0],[0,0,1],[-B[2, 1]*B[2, 2], -1, a1*B[2, 2]]]);

which avoids the Matrix product altogether.  But you can do better than that.  That is, you could also declare a module local alpha21 (the 2 means used in step 2), and assign it in Start

alpha21 := Matrix(3,3,[1,0,0],[0,0,1],[0,-1,0]]):

Then, in step 2, do

alpha21[3,1] := -B[2,1]*B[2,2];
alpha21[3,3] := a1*B[2,2];

Use alpha21 in step2 where you currently use alpha1.  The advantage of this is that it reuses the existing Matrix structure.  Doing this throughout your program will reduce the memory usage. 

Those work, but only because the row in a is originally 0.  He did mention adding to the row. 

Is there a reason you cannot use indexed names? Actually it isn't clear why you need names.  Could you just insert the results in a table (which is essentially what assigning them to an indexed name does?  Or an rtable (Array) if the indices are integers.

Why do you use the strange grouping in the number?  I write 10^8 as 100,000,000, not 1000,000,00.  Are you really calling fsolve that many times?   Even with fsolve(x=0) that will take about 8 hours on my machine.

You might consider using expressions rather than procedures.  That is,

y := x^2:
eval(y, x=2);
                       4
dy := diff(y,x);
                       2*x
eval(dy, x=2); 
                      4

There are trade-off with either technique, but it is easier to manipulate expressions.

Though not particulary useful, nor robust, here's a somewhat different approach.

a:=[["he","45",123,76,"1.0",4],["know","4",9,34,"3.2",5]];
parse(sprintf("%A",a));
                                [[he, 45, 123, 76, 1.0, 4], [know, 4, 9, 34, 3.2, 5]]

 V := Vector(convert([[1,6]], 'permlist', 6)): 
map(k -> V[k], {1,3,4});                     
                                                 {3, 4, 6}

You might want to look at the ?groups package

The catenation operator || does not evaluate its left operand.  Use cat.  However, the output of cat is not evaluated, so you'll need to do

eval(eval(cat(utility[3],1)), [alpha=1, beta=2]);

@corregumucio It still is not clear to me what you want to achieve.  However, you can manually simplify the expression. A nice way to see that is to freeze common subexpressions,

pw := piecewise(a < -(4*((-ce-pk)/b-(1/4)*(-ce-3*pk-2*ci)/b))*b
                , [{k <= (a-ce-pk)/b}]
                , And((4*((1/4)*(-ce-3*pk-2*ci)/b-(-ce-pk)/b))*b <= a, a < -(4*((-ce-pk)/b-(1/4)*(-ce-2*pk-2*ci)/b))*b)
                , [{k <= (1/4)*(3*a-ce-3*pk-2*ci)/b}]
                , (4*((1/4)*(-ce-2*pk-2*ci)/b-(-ce-pk)/b))*b <= a
                , [{k <= (1/4)*(3*a-ce-3*pk-2*ci)/b}]
                , []
               ):
# freeze all operands of relations.
pw1 := subsindets(simplify(pw), relation, rel -> map(freeze, rel));
            { [{k <= freeze/R1}]                  a < freeze/R4
            {
            { [{k <= freeze/R0}]        And(freeze/R4 <= a, a < freeze/R3)
     pw1 := {
            { [{k <= freeze/R0}]                  freeze/R3 <= a
            {
            {         []                            otherwise

From inspection it is clear that1 pw1 can be reduced to

pw2 := subsop(3=NULL,5=NULL,6=NULL,7=NULL,pw1);
                      { [{k <= freeze/R1}]        a < freeze/R4
               pw2 := {
                      { [{k <= freeze/R0}]          otherwise

Now thaw the frozen expressions

pw3 := thaw(pw2);
            {             a - ce - pk
            {      [{k <= -----------}]              a < 3 ce + pk - 2 ci
            {                  b
     pw3 := {
            {        3 a - ce - 3 pk - 2 ci
            { [{k <= ----------------------}]             otherwise
            {                 4 b

The following works

deqs := {NULL
         , diff(u2(y), y, y)+m*b*rho*h^2*GR*(c3*y+c4) = 0
         , diff(u1(y), y, y)+K*(diff(N(y), y))/(1+K)+GR*(c1*y+c2)/(1+K) = 0
         , -2*K*(2*N(y)+diff(u1(y), y))/(2+K)+diff(N(y), y, y) = 0
        }:

bcs := {NULL
        , u1(-1) = 0
        , u2(1) = 0
        , u1(0) = u2(0)
        , (D(u1))(0)+K*N(0)/(1+K) = (D(u2))(0)/(m*h*(1+K))
        , (D(N))(0) = 0
        , N(-1) = 0
       }:

dsolve(deqs union bcs);

The following works

plot(GAMMA, 1..2);
 
B := int(t^(a-1)*(1-t)^(b-1),t=0..1);    
                                                     GAMMA(b) GAMMA(a)
                                                B := -----------------
                                                       GAMMA(a + b)

plot3d(B, a=1..2, b=1..2);

Here is how to plot B vs a with b fixed

plot(eval(B, b=1.5), a = 1..2);

You can generate the tuples using

(r,l) := (4,2):                             
combinat:-permute([0$(r-l),1$r]);

If you are given T(x,t), then just compute a directly:

a := diff(T(x,t),t)/diff(T(x,t),x,x);

I see what you mean:

with(DifferentialGeometry): with(JetCalculus): with(Physics):
DGsetup([x], [u, psi], E, 1):
Setup(anticommutativeprefix={psi}):
f := psi[]*psi[1]:
TotalDiff(f,x);
                                 2
                           psi[1]  - psi[] psi[1, 1]

The problem lies in TotalDiff, it doesn't know anything about the anticommutation rule. That is, when it applies the Leibniz rule to differentiate the product, it uses standard multiplication rather than the overloaded multiplication operator. There is a way to hack around this, that is, you can replace the module local procedure used by TotalDiff with one that uses the overloaded multiplication operator.  An easy way to do that is to load the Physics package, so `*` is overloaded, convert the local procedure to a string, parse it, and then reassign it.  That works because the parsed multiply operation is now interpreted as the overloaded operation. Thus

kernelopts(opaquemodules=false): # needed to access the module local procedure
use td0=DifferentialGeometry:-JetCalculus:-TotalDiff:-td0 in 
    td0 := parse(sprintf("%a",eval(td0)));
end use:
TotalDiff(f,x);
                               -psi[] psi[1, 1]

No guarantees that this hack won't break something.

I'm not sure what the issue is.  Note that your linear transformation doesn't rotate the input, so the Maple plot won't look like the desired plot.  You can achieve that with a different transformation, say

A1 := <2,0;0,3>:
A2 := eval(<cos(theta),sin(theta);-sin(theta),cos(theta)>, theta=Pi/4):
A := A1 . A2;