restart;
#
# Define PDE
#
  pde:= E*M*diff(w(x,t),x$4)+A*rho*diff(w(x,t),t$2)=0;

#
# Assume a solution of the given form and substitute
# in the original pde to derive an ODE in 'x' only
#
  test1:=w(x,t)=W(x)*cos(omega*t+phi);
  pde1:=expand(eval(pde, test1)/(cos(omega*t)*E*M));
#
# Tidy up all the constants
#
  ode:=algsubs(A*rho*omega^2/(E*M)=alpha^4, pde1);
  sol:=simplify( convert(dsolve(ode), trig));
  sol:=subs( [ op([2,1,1],sol)=D1,
               op([2,2,1],sol)=D2,
               op([2,3,1],sol)=D3,
               op([2,4,1],sol)=D4
             ],
             sol
           );
#
# Try the first couple of boundary conditions
#
  b1:= simplify(rhs(eval( diff(sol, x), x=0))=0);
  b2:= simplify(rhs(eval( diff(sol, x$3), x=0))=0);
#
# Solve these BCs for S=D1, D3, and substitute the
# result back into the ODE solution
#
  bc1:=solve([b1, b2], [D1, D3])[];
  sol:=eval(sol, bc1);
 

#
# Now try the second pair of boundary conditions
#
   b3:= simplify(rhs(eval( diff(sol, x), x=L))/alpha=0);
   b4:= simplify(rhs(eval( diff(sol, x$3), x=L))/alpha^3=0);
#
# Consider
#
  b4-b3;
#
# This can only be true if either
#
#      D2=0, or
#      alpha*L=n*Pi (n an integer)
#
# Now consider
#
  b4+b3;
#
# This can only be true if either
#
#      D4=0, or
#      alpha*L=0
#
# I think(?) the only combination of the above two
# conditions which leads to non-trivial solutions
# is to have D4=0, and alpha*L=n*Pi. The latter
# condition however means restricting one of the
# variables in the definition of alpha - maybe
# omega??? Consider this later - for now, just make
# the substitutions D4=0, and alpha=n*Pi/L
#
  sol:=eval(sol, [D4=0, alpha=n*Pi/L]);
 

#
# Rebuild the complete solution
#
  ans:=eval(test1, sol);

#
# Notice that earlier, alpha=n*Pi/L was
# required, which in turn means that there
# is a restriction on omega, which can be
# determined as
#
  omega__n:= rhs~
             ( [ ( allvalues
                   ( isolate
                     ( A*rho*omega^2/(E*M)=(n*Pi/L)^4,
                       omega
                     )
                   )
                 )
               ]
             );
#
# Assuming negative frequencies are not allowed
# substitute the first of the omega__n values
# in the general solution
#
  ans:=eval(ans, omega=omega__n[1]);

  restart:
  ans := w(x) = _C1*exp(-I*(-K*E^3*i^3)^(1/4)*x/(E*i))+
                _C2*exp(-(-K*E^3*i^3)^(1/4)*x/(E*i))+
                _C3*exp(I*(-K*E^3*i^3)^(1/4)*x/(E*i))+
                _C4*exp((-K*E^3*i^3)^(1/4)*x/(E*i)):
  map2(op, 2, [op(rhs(ans))]);

restart;
L:= Vector( [ a*x^2+b*x*y-1,
              -(a*b-b)*x*y/a-z^2+(a-1)/a,
              -a*c*z^2/(b*(a-1))+(b+c)/b,
               a*x^2+b*x*y-1,
              -(a*b-b)*x*y/a-z^2+(a-1)/a,
              1
            ]
          ):
C:= Matrix( numelems(L),
            2,
            (i,j)-> `if`( j=1,
                          [ coeffs(L[i],[x,y,z], 't') ],
                          [t]
                        )
          );

  restart;
  PDE := (diff(u(x,y), x,x))-1/x*(diff(u(x, y), x)) -x^2*(diff(u(x,y), y,y))= 2;
  BCs := u(0, y) = 0, u(1, y) = 0:
  ICs := u(x, 0) = 0, D[2](u)(x, 0) = 2:
  num_sol := pdsolve(PDE, {BCs,ICs}, numeric);
#
# Access the procedure which will return the value
# of u(x,y) for any supplied values of the arguments
# x,y
#
  getUxy:=eval( u(x,y),
                num_sol:-value( u(x,y),
                                output=listprocedure
                              )
              );
#
# Now compute u(x,y) for any supplied x,y
#
  getUxy(0.1,2);
  getUxy(0.2,5);
  getUxy(0.3,7);
#
# Obviously you can compute/organise the data in a
# variety of ways - for example if you want a matrix
# of values for u(x,y) where rows correspond to
# x=0,0.1,0.2...0.5 and columns correspond to
# y=0,1,2..5, then you could use
#
  resMat:= Matrix
           ( 6,
             6,
             (i,j)-> getUxy((i-1)/10, j-1)
           );

  restart;
  ODE:=diff(y(x),x,x)+16*diff(y(x),x)=0:BCs:=y(0) =0, y(2)=3:
  ode_num_sol:=dsolve({ODE,BCs},numeric, output=listprocedure);
#
# Access the procedure which will return the value
# of y(x) for any supplied value of 'x'
#
  getY:=eval(y(x), ode_num_sol);
#
# Now compute the y(x) for any supplied x
#
  getY(0.1);
  getY(0.2);
  getY(0.5);
#
# Or if you wanted a matrix with x-values in column one, and
# associated y(x) values in column two, then you could use
#
  resMat:= Matrix
           ( 10,
             2,
             (i,j)->`if`( j=1,
                          (i-1.0)/10,
                          getY((i-1)/10)
                        )
           );

  restart;
#
# Alcuin's series. The coefficient of
# x^n ought to give the correct result.
# OP has n=36
#
  n:=36;
  expr:= x^3/((1-x^2)*(1-x^3)*(1-x^4)):
  ser:=taylor( expr, x=0, n+1):
  coeff(convert(ser, polynom), x, n);

#
# OP's calculation: code fine, but concept is wrong
#
  L := []:
  for a from 1 to 18 do
      for b from 1 to 18 do
          for c from 1 to 18 do
              if   `and`( a < b+c, b < a+c, c < b+a,
                          a > abs(b-c), b > abs(a-c), c > abs(a-b),
                          a+b+c = 36
                        )
              then L := [op(L), [a, b, c]]
              end if:
          end do:
      end do:
  end do;
  L;
  nops(L);

#
# Correct calculation. Note lower limits
# on for loops. Ensures that triples can
# only be obtained in increasing order, so
# avoiding the multiple counting
#
  L := []:
  for a from 1 to 18 do
      for b from a to 18 do
          for c from b to 18 do
              if   `and`( a < b+c, b < a+c, c < b+a,
                          a > abs(b-c), b > abs(a-c), c > abs(a-b),
                          a+b+c = 36
                        )
              then L := [op(L), [a, b, c]]
              end if:
          end do:
      end do:
  end do;
  L;
  nops(L);

plot( [ seq( Vc^(1/n)-0.5,
             n in [0.1,0.3,0.5,1.0, 3.0, 5.0,10.0]
            )
      ],
      Vc=0..1
    );
plots:-implicitplot( [ seq( Vc=(0.5+t)^n,
             n in [0.1,0.3,0.5,1.0, 3.0, 5.0,10.0]
            )
      ],
      Vc=0..1,
      t=-0.5..0.5
    );

plots:-implicitplot( [ seq( Vc=(0.5+t)^n,
             n in [0.1,0.3,0.5,1.0, 3.0, 5.0,10.0]
            )
      ],
      Vc=0..1,
      t=-0.5..0.5
    );