The bottleneck is not in Maple but rather in the access to the external site. Wrap your call in CodeTools:-Usage:

CodeTools:-Usage(X1()):
memory used=275.27KiB, alloc change=0 bytes, cpu time=0.00s, real time=0.59s

The cpu time (0.00) is the Maple kernel time.

While it won't help the time (it's hard to improve on 0), you can reduce the memory usage by using ?StringTools[StringBuffer] to catenate the lines:

X1 := proc ()
local str, sid, b, buf;
uses Sockets;
    buf := StringTools:-StringBuffer();
    sid := Open("google.com", 80);
    Write(sid, cat("GET /finance/historical?q=INDEXSP:.INX&histperiod=daily HTTP/1.0 \n\n"));
    do
        b := Read(sid);
        if b = false then break; end if;
        buf:-append(b);
    end do;
    Close(sid);
    buf:-value();
end proc:

CodeTools:-Usage(X1()):
memory used=38.13KiB, alloc change=0 bytes, cpu time=0.00s, real time=0.56s

I'm not familar with the details of the Maple plot command, however, I doubt that it is reentrant. Every evaluation of f, which is called by the main call to plot, attempts to create a new plot.

(**) f := z/(z-1):
(**) thaw(normal(subs(z=1/freeze(1/z), f)));
                                              1
                                         - --------
                                           -1 + 1/z

There are no differentials in that list of equations.  What are you attempting to solve for?

I don't know what you mean by "work on".  My immediate thought is to try what you suggested, that is, either assign q[O] as a piecewise expression,

q[O] := piecewise(x < 1, x, x^2):

or assign an equation to a variable:

eq := q[O] = piecewise(x < 1, x, x^2):

The latter may be more useful in that it allows you to use q[O] in other expressions.

I'll suggest saving the code as as a text file and using VerbatimInput from the fancyvrb LaTeX package.  This has the advantage in that any subsequent changes to the code do not require LaTeX edits.  Use the normal \section commands (after \appendix) in the LaTeX source. For example

\appendix
\section{Main Procedure}
\label{sec:maple-main}
\VerbatimInput{Main.mpl}
\section{Plot Routines}
\label{sec:maple-plot}
\VerbatimInput{Plot.mpl}

Consider using the separate text files as the source for the Maple package.  You can use the preprocessor $include directive in command-line Maple to handle stuff that is nested.  For example, your file layout may look something like

MyProg.mpl:

MyProg := module()
export ModuleLoad, Print, DoSomething;

$include <ModuleLoad.mpl>
$include <Print.mpl>
$include <DoSomething.mpl>

end module:
# save module to archive
LibraryTools:-Save(MyProg, "MyProg.mla"):

ModuleLoad.mpl:

ModuleLoad := proc()
...
end proc:

Running the following from the command-line creates the Maple archive MyProg.mla, which Maple automatically loads if it is located in the path defined by the ?libname variable.

$ maple MyProg.mpl

Did you look at your loop structure?  The inner loop is being executed more than 7777^5=28e18 times.

A few concepts have gone awry. The conditional statement, at the top of your post, checks the values of i and k before you have assigned them in the for-loop.  The loop itself keeps reassigning beta the same procedure.  The parameter i and k for beta have nothing to do with the loop variables (declared parameters to a procedure always have a local scope).  A simple fix is to eliminate the conditional (if) and the loops, and just assign alpha and beta the procedures you have.  In the call to Matrix, pass the values of i and k to alpha and beta.  That is

alpha := (i,k) -> ...
beta := (i) -> ...
Matrix(M, M, (i, k) →`if`(i = k, beta(i), alpha(i,k)));

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;