A symbolic solution is unlikely to exist, so you need to provide values
for your parameters to find the zeros numerically.

Here is a recipe how to do it for the real axis.

restart; Digits:=15;
 -(((beta*eta^2-(1/2)*beta+1)*p^2-(1/2)*beta^3)*KummerM((1/4)*((-beta+2)*
   p^2-beta^3)/p^2, 1, beta*eta^2)+     (1/2)*KummerM((1/4)*((-beta+6)*
   p^2-beta^3)/p^2, 1, beta*eta^2)*((beta-2)*p^2+beta^3))*
   exp(-(1/2)*beta*eta^2)/(p^2*beta);
 
 eval(%, [p = 1, eta = 1]); # some numerical values for parameters;
 
 f:=unapply(%, beta);

   beta -> 
     -((1/2*beta+1-1/2*beta^3)*KummerM(-1/4*beta^3-1/4*beta+1/2,1,beta)+
     1/2*KummerM(-1/4*beta^3-1/4*beta+3/2,1,beta)*(beta^3+beta-2))
     *exp(-1/2*beta)/beta

 # plot(f, -10 .. 10);

Maple does not see that there is a continuous extension to 0 being a zero.

 with(RootFinding);

  L:=NULL:
  t:=1e-3:               # starting point beyond beta=0
  for j from 1 to 3 do   # find 3 zeros, running to + infinity
    t:=NextZero(f, t);   # using that command, see the (crude) help
    L:= L, t;
  end do: t:='t':
  L:=[L];                # list of zeros

     L := [1.87588952396938, 2.58694685291078, 3.13796881735298]

 # check
 map(f, L);
                            -14                       -13
       [0.625993411295519 10   , -0.159057950051232 10   , -0.]

 # the negative values are zeros as well, a symmetry
 negL:=map(q -> -q, L);

  negL := [-1.87588952396938, -2.58694685291078, -3.13796881735298]

 # check that
 map(f, negL);
                        -14                      -13                       -13
  [-0.953311904713793 10   , 0.169099417324297 10   ,  0.336647449173820 10   ]

 # otherwise use g:=unapply(f(-x), x) and find zeros for g

Avoiding lists I would write that as

> convert(532, decimal, octal) -
> convert(235, decimal, octal);
> convert(%, octal);

                                 189

                                 275

 
> convert(275, decimal, octal) +
> convert(575, decimal, octal);
> convert(%, octal);

                                 570

                                 1072

> 2*8^0 + 7*8^1 + 0*8^2 + 1*8^3;

                                 570

f:=-144*z-44+12*sqrt(-12*z^3+96*z^2+24*z-15);
subs(z=x+I*y, f);
abs(%);
plot3d(%, x=-3..3,y=-3..3, axes=boxed);


Edited: 

Formally: if it has a root then bring the sqrt on one side and square it:

0=f;
% +144*z+44;
lhs(%)^2=rhs(%)^2;
lhs(%)-rhs(%);
expand(%);
#%/4; # gives Christian Wolinski's polynom
R:=solve(%);

       R := -4/3, -4/3, -4/3
       
The solution must be z = -4/3. But that is no solution:

eval(f, z=-4/3); simplify(%);

                                 296

NB: that is similar to ask for solving I = - sqrt(z).

2*ln(x+(x+ln(x))^(1/2))+2*(x+ln(x))^(1/2)

That will not work for several reasons, a better way is

H22*x*(cos(epsilon*((1/2)*L-x)/R)+kappa*cosh(epsilon*((1/2)*L-x)/R)):
int(%, x=0 .. L);

                        -30      2                      -12
    -2.95229265242607 10    omega  + 1.93138121006040 10  

You will lose that additional digits if you continue to operate with the result, but you can try something like this (without going into details for different length)

Digits:=10;
a:=.123456789;
m:=proc(a,b) evalf[2*Digits](a*b) end proc;

m(a,a);
                         0.015241578750190521

m(a,a) + 1;
                             1.015241579

I only have some test version, but there is a (general) workaround:

Go to the very directory, ..\bin.win,
mark the program file, maplew.exe, using the mouse 
and now use the context menu of the mouse (= right click) to add the program to the start menu.
Dito for the classical interface cwmaple.exe

I used a change x = 1/3000*46^(1/2)*z * t, t) assuming 0 < z and then
z = w / ( 1/3000*46^(1/2) ) to get

w^3*Int(t^2*ln(1+exp(-w*(t^2-1)^(1/2))),t = 0 .. infinity)

Now real / imaginary is decided in t=1 and using eval(Re(..)), 0 < z I arrive at

Int(1/2*t^2*ln((1+cos(w*(-t^2+1)^(1/2)))^2+sin(w*(-t^2+1)^(1/2))^2), t=0..1, 
  digits = 10, method = _d01ajc) +
Int(t^2*ln(1+exp(-w*(t^2-1)^(1/2))), t=1 .. infinity, 
  digits = 10, method = _d01amc);
w^3*%;
plot(%, w = 0 .. 30, numpoints=5);

That shows the rough locations of several extrema of a oscillating and
increasing function.

It does not directly answer your question, but for the above reformulation
you can diff for w formally and integrate and fsolve more "easily".

Hope that idea helps.

Here is mine, actually also handmade, but written in Maple (and perhaps
just some reformulation of the other answers):

  theAssumptions:= 0 < Omega, 0 < a, 0 < k, 0 < m;

Following Preben one can set Omega = 1 = m,

  subs(Omega=1, m=1, ex);
  P:=simplify(%) assuming theAssumptions;

    -2^(1/2)*(-a^2-2*k+a*(a^2+4*k)^(1/2))^(1/2) / ((a^2+4*k)^(1/2)-a)       

Maple has problems to find signs for the inputs to sqrt, but help it (and
have in mind strict convexity of sqrt):

  -a^2-2*k+a*(a^2+4*k)^(1/2):
  0 < (-1)*%;
  % +a*(a^2+4*k)^(1/2);
  lhs(%)^2 < rhs(%)^2; expand(%);
  is(%) assuming theAssumptions;

                         2             2       1/2
                    0 < a  + 2 k - a (a  + 4 k)

                                 true

Dito for the denominator. Now using sqrt(-p) = I*sqrt(p) for 0 < p we have

  P = -I *  2^(1/2)*(+a^2+2*k-a*(a^2+4*k)^(1/2))^(1/2)/((a^2+4*k)^(1/2)-a)

Maple still does not simplify the very term to 1, but knowing 'positive'
it now is enough to have it for the squares, like vv did.

  2^(1/2)*(+a^2+2*k-a*(a^2+4*k)^(1/2))^(1/2)/((a^2+4*k)^(1/2)-a);
  %^2; #map(expand, %);
  expand(numer(%)) / expand(denom(%)); # else denom is not expanded

                                  1

BTW: thx for the neat example

Often one way is to reduce to constant intervals.

In many cases order of integration does not matter (but you may check for your specific case), so interchange them and consider Int(Int(f(a,b),b=1-a..1+a),a=0..10). You may use Upper Case instead of lower case for integration untill the end.  

The inner integral can be mapped to any interval, in the following I map it to -1 ... 1. Then your problem has vanished.

  Int(f(a,b),b = 1-a .. 1+a);
  IntegrationTools:-Change(%, b = a*t+1, t): combine(%);
  Int(%,a=0..10);

                      10    1
                     /     /
                    |     |
                    |     |   f(a, a t + 1) a dt da
                    |     |
                   /     /
                     0     -1

Example:

  Int(Int(f(a,a*t+1)*a,t = -1 .. 1),a = 0 .. 10);
  eval(%, f = '(a,b) -> exp(-a^3-b)');
  value(%): subs(int=Int, %); # Maple does not know a symbolic solution
  evalf(%);

                   10    1
                  /     /
                 |     |         3
                 |     |   exp(-a  - a t - 1) a dt da
                 |     |
                /     /
                  0     -1

                10
               /
              |          3                  3
              |    exp(-a  + a - 1) - exp(-a  - a - 1) da
              |
             /
               0

                          0.370721316428104

Edited to add:

For that example good luck reduced the task even to dim = 1.

You can rewrite that final integrand as exp(-a^3)*exp(-1)*sinh(a)*2
to reduce (additive) cancellation errors if you want high precision
(though Maple should care for that on its own).

It is like integrating in a pole:

Int((1/2)/(s^(1/2)*GAMMA(1/2)*(t^(1/2)-s^(1/2)))*(s^6), s = 0 .. t);
IntegrationTools:-Change(%, s=x*t, x);
IntegrationTools:-GetIntegrand(%);
expand(%);
series(%, x=1, 3);

       2/(x-1)+23/2+219/8*(x-1)+O((x-1)^2)

So it is like int(1/(x-1), x=0..1), since the numerator in x=1 is finite.

You want to expand it (the ordering may change)

  x= (a+b*c)/d: 
  expand(%);
                                b c
                            x = --- + a/d
                                 d

It seems that 'simplify' writes (u-1)*(u+1) as u^2 - 1 and then does not factor again

restart: interface(version); # mine is Maple 18
c:= (1+sqrt(5))/2: #c:=2;
theAssumption:= 1<u,u<c:

g:=((u-1)*(u+1))^(1/2)*u/(u*(u-1))^(1/2);

simplify(g);

                              2     1/2
                            (u  - 1)    u
                            --------------
                                       1/2
                            (u (u - 1))

combine(g) assuming theAssumption;
simplify(%) assuming theAssumption;
#combine(%) assuming theAssumption;

                                  1/2  1/2
                           (u + 1)    u

In your case you can scale the matrix by it's determinant, solve and scale back. Find a file attached

MP_EigenValues_scaling.mw

MP_EigenValues_scaling.pdf

numer(%)/combine(denom(%)) and use Pi*x=arccos(r), Pi*y=arccos(s),
use cos(n*arccos(r)) = orthopoly[T](n,r) to come up with

Int(Int(
  2/(-s^2+1)^(1/2)/Pi^2*(T(n,r)*T(n,s)-1)/(-2+r+s)/(-r^2+1)^(1/2),
  r = -1 .. 1),s = -1 .. 1);

         1    1
        /    /
       |    |            2 (T(n, r) T(n, s) - 1)
       |    |   ------------------------------------------ dr ds
       |    |      2     1/2   2                 2     1/2
      /    /    (-s  + 1)    Pi  (-2 + r + s) (-r  + 1)
        -1   -1

Do it stepwise to get results for numerical values of n.