The problem is that once you refer to parameter by name in the procedure, the corresponding argument must have been passed, otherwise an error is raised. One solution is to query nargs to determine whether change was passed.  That is

   if nargs = 0 then currentdir("f:/") else ... end if

A somewhat neater solution is to give change a default value

changedir := proc( change :: string := NULL, $)
    if change = NULL then ... end if
end proc:

Note that it is usually legal to assign a default value (here NULL) that is different from the declared type of the parameter. 

Why?  There may be a better way to do what you want.  I assume you mean a Maple ?table, and not an ?rtable, which includes Matrices and Arrays).  There isn't a really good way to do so, that is, one that is efficient.  One way is nops([indices(mytable)]). The reason that is not ideal is that the calll to ?indices creates a sequence of all the indices in the table; that uses memory.  Possibly not a concern here. 

The I is the imaginary unit, that is, sqrt(-1).

Is an example of what you want?

dsys := { NULL
          , diff(x1(t),t,t) = -x1(t)
          , diff(x2(t),t,t) = -x2(t)
          , x1(0) = 0, D(x1)(0)=1
          , x2(0) = 1, D(x2)(0)=0
        }:
integ := dsolve(dsys, numeric, 'events' = [[x1(t)=x2(t) , [x1(t)=x1(t)+x2(t), x2(t)=1]]]):
plots:-odeplot(integ, [[t,x1(t)],[t,x2(t)]], 0..3*Pi);

Note that while you can change the values of the state variables (x1,x2) at an event, you cannot change the form of the differential equations. You could, however, change a discrete-time variable that causes a branch in a piecewise expression to change, which effectively changes the form of the equations.  For example,

dsys := { NULL
          , diff(x1(t),t) = piecewise(b(t)=0, 2-x1(t), -2+x1(t)^2)
          , x1(0) = 0
          , b(0) = 0
        }:

integ := dsolve(dsys, numeric
                , 'events' = [[x1(t)=1, b(t)=1]]
                , 'discrete_variables' = [ b(t) :: boolean ]
               ):
plots:-odeplot(integ, [[t,x1(t)]], 0..2);

Using differential forms would make life a bit simpler.  That is, what you are really being asked (with forms being disguised as vector fields) is, given a two-form B, find a one-form A such that dA = B, with dA the exterior derivative of A. Clearly, this notation suggests an integration is in order. The additional constraint you want is dB=0. Because B is exact, that is, it is the exterior derivative of a one-form, it must be closed, that is, its exterior derivative must be zero.

The ?DifferentialGeometry package can handle these computations, however, learning to use it may be daunting. Here it is, regardless.  Possibly it will inspire you to learn about differential forms; the notation is clean and they work in all dimensions (unlike the vector calculus, which is restricted to to three).

with(DifferentialGeometry):
DGsetup([x,y,z],M):

B := evalDG(Bx(x,y,z)*dy &w dz + By(x,y,z)*dz &w dx + Bz(x,y,z)* dx &w dy):

# Here are the integrability conditions
dB := ExteriorDerivative(B);
        //d             \   /d             \   /d             \\
  dB := ||-- Bx(x, y, z)| + |-- By(x, y, z)| + |-- Bz(x, y, z)|| dx ^ dy ^ dz
        \\dx            /   \dy            /   \dz            //

eq := Tools:-DGinfo(dB, "CoefficientSet")[] = 0;
             /d             \   /d             \   /d             \
       eq := |-- Bx(x, y, z)| + |-- By(x, y, z)| + |-- Bz(x, y, z)| = 0
             \dx            /   \dy            /   \dz            /

# To find A such that dA = B, we essentially integrate B along paths.
# The DeRhamHomotopy procedure does this:
DeRhamHomotopy(B);
   1
  /
 |
 |   _z1 (z By(_z1 x, _z1 y, _z1 z) - y Bz(_z1 x, _z1 y, _z1 z)) d_z1 dx +
 |
/
  0

       1
      /
     |
     |   _z1 (-z Bx(_z1 x, _z1 y, _z1 z) + x Bz(_z1 x, _z1 y, _z1 z)) d_z1 dy
     |
    /
      0

     -

       1
      /
     |
     |   _z1 (-y Bx(_z1 x, _z1 y, _z1 z) + x By(_z1 x, _z1 y, _z1 z)) d_z1 dz
     |
    /
      0


# Now take an explicit example
B := evalDG(x^2*dy &w dz + (y+x)*dz &w dx -(2*x+1)*z* dx &w dy);
                                                           2
           B := -(2 x + 1) z dx ^ dy + (-y - x) dx ^ dz + x  dy ^ dz


# Verify that its exterior derivative is 0
ExteriorDerivative(B);
                                 0 dx ^ dy ^ dz


# Compute A
A := DeRhamHomotopy(B);
                                                   2
A := (1/2 z y x + 2/3 z y + 1/3 z x) dx + (-3/4 z x  - 1/3 z x) dy

               2                  2
     + (1/4 y x  - 1/3 x y - 1/3 x ) dz


# Verify B = dA
dA := ExteriorDerivative(A);
                                                            2
           dA := (-2 z x - z) dx ^ dy + (-y - x) dx ^ dz + x  dy ^ dz

Tools:-DGequal(B, dA);
                                      true

You can use ?frem for floating points.  See the help page for details.

An Array might be a decent choice:

rus := Array(1..3, [0,0.2,0.5]);
for j from 15 to 27 by 3 do
    for k to 3 do
           ru:= rus[k];
    end do
end do;

You could also use a list

rus := [0, 0.2, 0.5]:
for j from 15 to 27 do
    for ru in rus do
    end do;
end do;

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