1. Each of your loops runs from 1 to 6 but the functions U[index], V[index], W[index] are only defined for index =1..2
2. to perform integrations numerically, it is much more efficient to use evalf(Int(Int(...))), rather than evalf(int(int(...))). Not the capitalisation of the command Int()

The attached generates values for for the 8 integrals which are defined for N=2, reasonably quickly (a few secs on my machine)

someInts.mw

provided that you only want something which "looks pretty"

restart;
for i to 2 do
    for j to 2 do 
        S[i, j] := 2*mu*varepsilon[i, j] - add( (2*mu)/3*varepsilon[r, r]*delta[i, j],
                                                r = 1 .. 3
                                              );
        print
        ( parse
          ( cat("S[", i, j, "]")
          )
          =
          parse
          ( StringTools:-Subs
            ( [ "," = "" ],
              convert
              ( S[i,j],
                string
              )
            )
          )
        ); 
     end do;
end do;

which produces

is shown  by

  restart:
  R:=10:
  r:=5:
  plot3d( [ ( R+r(u,v)*cos(v) )*sin(u),
            ( R+r(u,v)*cos(v) )*cos(u),
            r*sin(v)
          ],
          u=0..2*Pi,
          v=0..2*Pi,
          style=patchnogrid,
          scaling=constrained,
          coordinateview=[-15..15, -15..15,-5..5],
          axis[3]=[tickmarks=[3, subticks=4]],
          lightmodel=light3,
          viewpoint = circleleft
        );

which produces

are shown in the code below

  restart;
  h:= 0.03;
  C:= 1;
  rho:= 7800;
  omega[m, n] := 222;
  S:=(Matrix([[rho*h*omega[m,n]^2+100,200],[280,rho*h*omega[m,n]^2+444]])).(Vector(1..2,[[A],[B]]))=-C*(Vector(1..2,[555,6665]));
#
# One way
#
  solve([entries(lhs(S)=~rhs(S),nolist)]);
#
# Another way - setting up the problem slightly
# differently
#
  restart;
  with(LinearAlgebra):
  h:= 0.03;
  C:= 1;
  rho:= 7800:
  omega[m, n]:= 222:
  M:=Matrix([[rho*h*omega[m,n]^2+100,200],[280,rho*h*omega[m,n]^2+444]]);
  V:=-C*Vector(1..2,[555,6665]);
  solve(GenerateEquations(M, [A,B], V));

see the attached

sol.mw

to the question here

https://www.mapleprimes.com/questions/231383-Hexagon-Inscribed-And-Exinscrit-In-A-Circle,

then the code

  restart;
  with(geometry):
  with(plots):
  _EnvHorizontalName:= x:
  _EnvVerticalName:= y:
  local gamma:
  R:=3.0:
  angs:=Array( 0..5,
               [0, (1/3)*Pi, 3*Pi*(1/4)+1/5, 7*Pi*(1/6)+2/5, 8*Pi*(1/5), 13*Pi*(1/7)]
             ):
#
# Define the origin OO
#
  point(OO, [0, 0]):
#
# Define the circle
#
  circle(C, [OO, R]):
#
# Define points on the circle as specified by the
# entries in the array 'angs'
#
  seq( point( cat(P, j), [R*cos(angs[j]), R*sin(angs[j])]), j=0..5):
#
# Generate the tangents at these points
#
  seq( TangentLine( cat(T, j), cat(P, j), C),
       j=0..5
     ):
#
# Get the points of intersection of the tangents
#
  seq( intersection( cat(Q, j), cat(T, j), cat(T, irem(j+1, 6) ) ),
       j=0..5
     ):
#
# Produce the line segments which form the hexagon
#
  seq( segment( cat(S, j), cat(Q, j), cat(Q, irem(j+1,6))),
       j=0..5
     ):
#
# Produce the lines through adjacent tangent points
#
  seq( line( cat(L, j), [cat(P, j), cat(P, irem(j+1,6))]),
       j=0..5
     ):
#
# Produce the lines through the hexagon vertices
#
  seq( line(cat(V, j), [cat( Q, j), cat(Q, j+3)]), 
       j=0..2
     ):
#
# The OP's points alpha, beta, gamma - not sure what to
# call these!
#
  pts:=[beta, gamma, alpha]:
  seq( intersection( pts[j+1], cat(L, j), cat(L, j+3)), j=0..2):
  line(M, [alpha, gamma]):
#
# Draw the geometric structures
#
  p1:= draw( [ C(color=black),
               seq( cat(S, j)(color=green, thickness=3), j=0..5),
               seq( cat(L, j)(color=magenta, linestyle=spacedash), j=0..5),
               seq( cat(V, j)(color=red), j=0..2),
               M(color=blue, thickness=3)
             ],
             view=[-6..10, -15..6],
             axes=none,
             scaling=constrained,
             size=[800, 800]
           ):
#
# Place text at appropriate locations
#
  p2:= textplot( [ [ coordinates(alpha)[], "alpha", align={above, left}, font=[times, bold,18]],
                   [ coordinates(beta)[],  "beta",  align={above, left}, font=[times, bold, 18]],
                   [ coordinates(gamma)[], "gamma", align={above, left}, font=[times, bold, 18]],
                    seq( [coordinates(cat(P,i))[], cat("P", i),align = [above, left] ], i = 0 .. 5),
                    seq( [coordinates(cat(Q,i))[], cat("Q", i),align = [above, left] ], i = 0 .. 5)
                 ],
                 axes=none
               ):
#
# Display geometry and text
#
  display( p1,p2);

will produce the output

I didn't do this in my original response, because I thought it made the plot a bit too "cluttered"

See the attached

hex2.mw

MAny of the calculations which you are performing can be performed more siimply by using commands from the geometry package.

The code

  restart;
  with(geometry):
  with(plots):
  _EnvHorizontalName:= x:
  _EnvVerticalName:= y:
  local gamma:
  R:=3.0:
  angs:=Array( 0..5,
               [0, (1/3)*Pi, 3*Pi*(1/4)+1/5, 7*Pi*(1/6)+2/5, 8*Pi*(1/5), 13*Pi*(1/7)]
             ):
#
# Define the origin OO
#
  point(OO, [0, 0]):
#
# Define the circle
#
  circle(C, [OO, R]):
#
# Define points on the circle as specified by the
# entries in the array 'angs'
#
  seq( point( cat(P, j), [R*cos(angs[j]), R*sin(angs[j])]), j=0..5):
#
# Generate the tangents at these points
#
  seq( TangentLine( cat(T, j), cat(P, j), C),
       j=0..5
     ):
#
# Get the points of intersection of the tangents
#
  seq( intersection( cat(Q, j), cat(T, j), cat(T, irem(j+1, 6) ) ),
       j=0..5
     ):
#
# Produce the line segments which form the hexagon
#
  seq( segment( cat(S, j), cat(Q, j), cat(Q, irem(j+1,6))),
       j=0..5
     ):
#
# Produce the lines through adjacent tangent points
#
  seq( line( cat(L, j), [cat(P, j), cat(P, irem(j+1,6))]),
       j=0..5
     ):
#
# Produce the lines through the hexagon vertices
#
  seq( line(cat(V, j), [cat( Q, j), cat(Q, j+3)]), 
       j=0..2
     ):
#
# The OP's points alpha, beta, gamma - not sure what to
# call these!
#
  pts:=[beta, gamma, alpha]:
  seq( intersection( pts[j+1], cat(L, j), cat(L, j+3)), j=0..2):
  line(M, [alpha, gamma]):
#
# Draw the geometric structures
#
  p1:= draw( [ C(color=black),
               seq( cat(S, j)(color=green, thickness=3), j=0..5),
               seq( cat(L, j)(color=magenta, linestyle=spacedash), j=0..5),
               seq( cat(V, j)(color=red), j=0..2),
               M(color=blue, thickness=3)
             ],
             view=[-6..10, -15..6],
             axes=none,
             scaling=constrained,
             size=[800, 800]
           ):
#
# Place text at appropriate locations
#
  p2:= textplot( [ [ coordinates(alpha)[], "alpha", align={above, left}, font=[times, bold,18]],
                   [ coordinates(beta)[],  "beta",  align={above, left}, font=[times, bold, 18]],
                   [ coordinates(gamma)[], "gamma", align={above, left}, font=[times, bold, 18]]
                 ],
                 axes=none
               ):
#
# Display geometry and text
#
  display( p1,p2);

will produce the same plot as your worksheet.

See the attached worksheet

hex.mw

I can reproduce this by setting

Tools->Options->Display->Ouput Display to "Maple Notation"

In general this entry should be "2-D Math Notation"

and try to define your system in some kind of vaguely consistent way. Your posted code states

pde := W*diff(C(z, t), z $ 2) - diff(C(z, t), t $ 1) = 0
bc[1] := C(x, 0) = 0;
bc[2] := D[1](f)(0, t) = VMK*D[1](f)(t)/ZRTWA;
bc[3] := C(infinity, t) = 0;

In your first boundary condition, ie

bc[1] := C(x, 0) = 0;

what is the name 'x'? Where else in this problem does it appear? Is it just some other (random) independent variable? Or maybe you meant 'z'?

In your second boundary condition, ie

bc[2] := D[1](f)(0, t) = VMK*D[1](f)(t)/ZRTWA;

Why are you setting a boundary condsion on a function 'f', when no function 'f' occurs anywhere in your pde system? Is this a boundary condition on the function 'C?

Now I could probably guess at what you actually meant to type in the above two situations - but frankly I'm having a bad day and I see no reason why I should have to make educated guesses, just because you can't be bothered to define your problem accurately

for me to fix with confidence - In places, I'm making wild guesses about what you intended.

The attached now "works" and produces answers, but you are going to have to check it line by line to ensure that it is what you want.

A non-exhaustive list of things I fixed

1. You use functions n() and c(), together with parameters 'n' and 'c' - don't. IIRC I changed all instances of the parameter 'n' to 'nn' and 'c' to 'cc'
2. You need 6 bcs. You have 6 but two of them just repeat conditions you already have: in other words you actually have four - so I made up a couple of "new" bcs, just for testing purposes
3. Your loop relies on varying a quantity 'bval' which does not occur anywhere in your pdes, bcs - so even if your code functioned, nothing depends on 'bval' and so you would get the same answer every time. I assumed you were trying to vary the parameter beta1, so I fixed the definition of parameters so that this happens correctly.
4. Almost certainly several other things which I have now forrgotten

The attached now runs, and produces plots, although the dependence on the value of the parameter 'beta1' is very small - so small that it is almost impossible to detect using the sol[n]:-plot3d() method, but it is (just) detectable using the  sol[n]:-value() method.

Please check the attached very carefully - because I have changed so much, it is almost inevitable that I got something wrong somewhere. In particular verify the bcs I added

sim.mw

use algsubs() rather than subs()

See the attached

Download subsProb.mw

Don't worry too much about the upload error on this site. Someone at Maple broke this about a month ago and depite numerous complaints still can't be bothered to fix it!!

there are numerous typos(?) which probably mean that the integral(s) will never be calculated!

If I make guesses as to what you intended (every single one of which you should check, because I am guessing)  then the attached will run code will run and produce sensible(?) answers.

Note that I have used N=3 in the attached - so 3^4=81 integrals are being performed and this takes about 40secs on my machine

With N=6, 6^4=1296 integrals will be performed, so on my machine this would (probably) take around 640secs (call it 10minutes) on my machine - I can't be bothered to wait that long - maybe you can!!

oddInt.mw

You read the error message which Maple provides telling you that you have a recursives assignment statement.

You then eliminate this recursive assignment statement, since it can never be resolved

The offending stetement is highlighted in the attached

restart:

Define the parameters
N:=0.1:
f := l -> 1/(a*sqrt(2*Pi))*exp(-l^2/(2*a^2)):
alpha:=0.1:
phi:=0:
M:=sqrt(N*(N+1))*exp(I*phi):
g:=0.1:
Gamma:=g/2*(1+2*N):
delta:=0.25:
Delta:=0:
k:=alpha+I*delta:
Omega:=0.2:
z:=-0.5;
x:=-0.0120296652120473-0.389736176167980*I:
a:=2.5:
sum1[-10]:=0:
h:=0.1:

                           z := -0.5


for n from -10 to 10 do
    d[n]:=n*h;
    w[n]:=n*h:
    s[n]:=I*d[n]:
    q[n]:=delta+Delta-w[n]:
    if   N*(N+1)>(delta+Delta-w[n])^2
    then r := w -> sqrt(N*(N+1)-(delta+Delta-w)^2):
    elif N*(N+1)<(delta+Delta-w[n])^2
    then r := w -> I*sqrt(-(N*(N+1)+(delta+Delta-w)^2)):
         a1[n]:=(r(w[n])+I*q[n])/(s[n]+Gamma-r(w[n])+I*delta)+(r(w[n])-I*q[n])/(s[n]+Gamma+r(w[n])+I*delta);
         a2[n]:=((r(w[n])+I*q[n]+2*alpha)/(alpha*(alpha+r(w[n])))+g*M/(conjugate(k)*(r(w[n])+conjugate(k)))):
         k2[n]:=(r(w[n])+I*q[n])/(Gamma+r(w[n])+conjugate(k)):
         k1[n]:=-1/(4*r(w[n])*(s[n]+Gamma-r(w[n])+I*delta)):
         k3[n]:=g*conjugate(M)/(Gamma+r(w[n])+k)*((r(w[n])+I*q[n]+2*k)/(k*(r(w[n])+k))+g*M/(alpha*(alpha+r(w[n])))):
         c1[n]:=1/(4*r(w[n])*(s[n]+Gamma+r(w[n])+I*delta)):
         c2[n]:=(r(w[n])-I*q[n])/(Gamma-r(w[n])+conjugate(k)):
         b2[n]:=((r(w[n])-I*q[n]-2*alpha)/(alpha*(alpha-r(w[n])))+g*M/(conjugate(k)*(r(w[n])-conjugate(k)))):
         c3[n]:=g*conjugate(M)/(Gamma-r(w[n])+k)*((r(w[n])-I*q[n]-2*k)/(k*(r(w[n])-k))+g*M/(alpha*(alpha-r(w[n])))):
         f1[n]:=1/(4*r(w[n])*alpha*(r(w[n])-alpha)*(s[n]+Gamma-r(w[n])+2*alpha+I*delta))*((r(w[n])+I*q[n])*(r(w[n])+I*q[n]-2*alpha)/(Gamma+r(w[n])-k)-g^2*abs(M)^2/(Gamma+r(w[n])-conjugate(k)))+1/(4*r(w[n])*alpha*(r(w[n])+alpha)*(s[n]+Gamma+r(w[n])+2*alpha+I*delta))*(
(r(w[n])-I*q[n])*(r(w[n])-I*q[n]+2*alpha)/(Gamma-r(w[n])-k)+g^2*abs(M)^2/(Gamma-r(w[n])-conjugate(k)))+2*((k-Gamma+I*q[n])/((Gamma-k)^2-r(w[n])^2)+(g*M)/((Gamma-conjugate(k))^2-r(w[n])^2))+(Gamma+conjugate(k)-I*q[n])/(((Gamma+conjugate(k))^2-r(w[n])^2)*(s[n]+2*Gamma+conjugate(k)+I*delta))-g*conjugate(M)/(((Gamma+k)^2-r(w[n])^2)*(s[n]+2*Gamma+k+I*delta))+g*M/(4*r(w[n])*conjugate(k))*((r(w[n])+I*q[n]-2*conjugate(k))/((Gamma+r(w[n])-conjugate(k))*((r(w[n])-conjugate(k))*(s[n]+Gamma+2*conjugate(k)-r(w[n])+I*delta)))-(r(w[n])-I*q[n]+2*conjugate(k))/((Gamma-r(w[n])-conjugate(k))*((r(w[n])+conjugate(k))*(s[n]+Gamma+2*conjugate(k)+r(w[n])+I*delta))))+g*conjugate(M)/(4*r(w[n])*k)*((r(w[n])+I*q[n])/((Gamma+r(w[n])-k)*((r(w[n])-k)*(s[n]+Gamma+2*k-r(w[n])+I*delta)))-(r(w[n])-I*q[n])/((Gamma-r(w[n])-k)*((r(w[n])+k)*(s[n]+Gamma+2*k+r(w[n])+I*delta)))):
         cp[n]:=1/(2*r(w[n]))*a1[n]+Omega^2*((k1[n]*(k2[n]*a2[n]+k3[n])+c1[n]*(c2[n]*b2[n]+c3[n]))+f1[n]):
         cz[n]:=-2*I*Omega*(1/(2*r(w[n])*(s[n]+Gamma+r(w[n])+I*delta))*((r(w[n])-I*q[n])/(Gamma+conjugate(k)-r(w[n]))-g*conjugate(M)/(Gamma+k-r(w[n])))+1/(2*r(w[n])*(s[n]+Gamma-r(w[n])+I*delta))*((r(w[n])+I*q[n])/(Gamma+conjugate(k)+r(w[n]))+g*conjugate(M)/(Gamma+k+r(w[n])))-(Gamma+conjugate(k)-I*q[n])/(((Gamma+conjugate(k))^2-r(w[n])^2)*(s[n]+2*Gamma+conjugate(k)+I*delta))+g*conjugate(M)/(((Gamma+k)^2-r(w[n])^2)*(s[n]+2*Gamma+k+I*delta)));             
         L[n]:=(-2*cp[n]*z+cz[n]*x)*f(w[n]);
#
# The following is a recursive assignment. The first
# time the loop code executes, (ie n=-10), then this
# statement is OK because you have defined sum1[-10]
# before entering the loop.
#
# For the next loop index value you are trying to
# assign sum1[-9] in terms of sum1[-9], which is undefined
# So the assignment is recursive
#
          sum1[n]:=(L[n]+sum1[n])/3;

end if:
od:

The most important thing is to pick the correct coordinate system - and understand it.

The code below is just one of many possibilities

plot3d(3+cos(h), theta=0..2*Pi, h=-1..1, coords=cylindrical, style=surface, axes=none);
plot3d(3-cos(h), theta=0..2*Pi, h=-1..1, coords=cylindrical, style=surface, axes=none);

produces

and

the issue you are having.

Consider the code

#
# Create a 1D Array each of whose elements is a
# a list containing a sublist and a value, ie are
# of the form
#        [[x,y], z]
#
   A1:=Array(1..10, [seq([[i,i+1],i^2],i=1..10)]):
#
# Check an element
#
  A1[3];
#
# Extract 'pieces' from this element
#
  A1[3][1];
  A1[3][2];
#
#
# Create a 2D Array, where entries in the first
# "column", ie A2[n,1] are lists and entries in
# second "column" are simple values
#
  A2:=Array(1..10, 1..2,[seq([[i,i+1],i^2],i=1..10)]):
  A2[3,1];
  A2[3,2];

which will produce the output

                          [[3, 4], 9]

                             [3, 4]

                               9

                             [3, 4]

                               9

How you choose to arrange data within an array and the number of dimensions you might want, depend almost entirely on your application and how you subsequently wish to access the data.

The attached worksheet shows

1. For your first example, both analytic and numeric solutions are available. The relevant execution group shows
   1. A plot of u(x, 1/2) with x=0..1, for both analytic and numeric solutions. (In the analytic solution, only the first 100 terms of the infinite sum are used)
   2. A plot of u(1/2, t) with t=0..1, for both analytic and numeric solutions. (In the analytic solution, only the first 100 terms of the infinite sum are used)
   3. A plot of u(x ,t) with x=0..1,  t=0..1,for both analytic and numeric solutions. (In the analytic solution, only the first 100 terms of the infinite sum are used)
2. For your second example, no analytic solution is available. To obtain a numeric solution it is necessary to assign values for the parameters a, f,: k, v, L, T__0, Abitrarily I set all of these equal to 1. (Feel free to use other values!) The worksheet shows
   1. A plot of u(x, 1/2) with x=0..1, for the numeric solution
   2. A plot of u(1/2,  t) with t=0..1, for the numeric solution
   3. A plot of u(x ,t) with x=0..1,  t=0..1,for the numeric solution

pdePlots.mw