I'm willing to explain what I can, but Maple doesn't currently have SOC programming.  Do you want to write your own program for this?  It seems like a nontrivial task.  Or do you just want to use the Optimization package as-is?
This throws away the efficiency benefits of SOC programming.

That's strange.  I didn't get any lost connection with the number of loops set to 10.  The results were

                            266.6730
                            248.3509
                            236.7313
                            218.5128
                            148.3256
                            277.3509
                            234.7313
                            243.5128
                            74.3256
                            203.3509

This is using Maple 14.01 under Windows 7.  What Maple version and operating system are you using?

Look up "digital signature", w.g. on Wikipedia

What initial conditions are supposed to produce those curves? 

What do you mean by "proper"?  Are you saying that those are the actual differential equations, or are you referring to some mathematical property of the equations?  What theory are you referring to?  How were those graphs produced?

It looks to me like RegSubs may not be using a very efficient algorithm.  But as I understand it you are replacing a single fixed pattern "\n     #" by "&\n      ", so you should be able to use SubstituteAll instead of RegSubs.

You didn't show what the curves should look like.  The differential equations you're solving are

-diff(beta(t),t)+1.428571428*sin(beta(t))*sin(teta(t))/(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t)))+.2500000000*tan(beta(t))*cos(beta(t)+ang(t))/(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))) = 0,

-diff(ang(t),t)-1.428571428*sin(beta(t))*cos(teta(t))/(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t)))+.2500000000*tan(beta(t))*sin(beta(t)+ang(t))/(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))) = 0,

-diff(teta(t),t)-1.428571428*sin(beta(t))*(sin(teta(t))*sin(beta(t))+.35*sin(beta(t)+ang(t))*sin(teta(t))+cos(teta(t))*cos(beta(t))+.35*cos(teta(t))*cos(beta(t)+ang(t)))/(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t)))-.2500000000*tan(beta(t))*(cos(beta(t)+ang(t))*sin(beta(t))-sin(beta(t)+ang(t))*cos(beta(t)))/(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t)))+.5*cos(beta(t)) = 0

Are these correct?  With the initial conditions beta(0) = 0, ang(0) = 0, teta(0) = 0 the solution is

beta(t)  = 0, ang(t) = 0, teta(t) = t/2

and this is what your numerical solutions produce.  With beta(0) = 0.01, the solution appears to have a singularity near t = 4.7655932, as the denominator cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t)) of the left sides of your differential equations goes to 0.  If you used the default numerical methods, Maple would tell you this with

Warning, cannot evaluate the solution further right of 4.7655932, probably a singularity

As far as I can tell, Maple is correct: there really is a singularity in the solutions near that point.  If that does not fit with what you expect, the error is likely to be in your equations.

Using method=classical[foreuler] avoids the warning, but the solution it produces for t > 4.7655932 has no connection with reality.

@amrramadaneg : I think you want k = i+j-1 mod M + 1 and l = i-j-1 mod M + 1..

"Left precedence" refers to things that are on the left of mod, and "right precedence" refers to things that are on the right of mod.  So in a + b mod c + d, the + on the left has higher precedence than the mod, and the + on the right has lower precedence than the mod.  Thus the result is ((a+b) mod c) + d.

When you ask to solve for just d, you're asking for a value of d that will make both equations true (as functions of the parameters r, rr, a, p).  There is none, so Maple returns nothing.

Normally with two equations, you would expect to solve for two variables.  So e.g.

> solve({r*rr=0, a*p*d=rr},{d, rr} );

   {d = 0, rr = 0}

T does not contain x, because you substituted (1/2)*a for x.  It doesn't have any terms containing U(t) either.  It is just
-diff(U[2,1](t),t)^2*sinh(pi)^2*sin(pi*y/b)^2.  By the way, you probably want to use Pi rather than pi: Pi is the constant 3.14159..., pi is just another variable.  As for your question, if you want to differentiate with respect to diff(U[2,1](t),t), the simplest way is to temporarily substitute a new variable:

> Tv:= subs(diff(U[2,1](t),t) = v, T);
   diff(Tv,v);

and then if you wish you can substitute back:

> subs(v = diff(U[2,1](t),t), %);

Somewhat simpler perhaps:

M:= 100:
A:= Array(1..2, 1..15, 1..M, (s,i,j) -> i + (-1)^s * j  mod M);

I could be wrong, but I doubt that the inverse laplace transform of exp(x*sqrt((a*s^2+b*s)/(1+c*s^2))) can be expressed in closed form.  Do you have any reason to think that it can?

Yes, or alternatively JacobiTheta3(z,q) = sum(q^(n^2)*exp(2*I*n*z), n=-infinity .. infinity)

and thus for real nonzero sigma,

The sum for n = -infinity to infinity can be expressed using a Jacobi theta function, but for some reason Maple does not know about such sums.

.

Ah, I see: orthopoly[P](n,x) is an ordinary polynomial in x, and the coefficients get large, causing numerical difficulties. I guess the use of LegendreP allows more stable methods of computation.