Solving for complex phase speed eigenvalue

[primalfido]

Hey there, first timer poster here

Homework Statement

I'm working on a barotropic linear instability analysis and I've been having trouble getting an expression for the complex phase speed eigenvalue $C = C_r + i*C_i$ for the purpose of plotting a dispersion diagram ($C_i$ vs $k$ or $C_r$ vs $k$)

Homework Equations

The Normal mode solution:

$\psi(x,y,t) = \Psi(y)e^{ik(x - Ct)}$

where k is the zonal wavenumber

And the eigenvalue problem I'm trying to solve for C

$\frac{\partial^{2} \Psi}{\partial y^{2}} \Psi + (\frac{\beta - \frac{\partial^{2} U}{\partial y^2}}{U-C} -k^{2} )\Psi = 0$

where U is a velocity profile that's given in gaussian form \e^{-y^2} and \beta is a given constant.

The Attempt at a Solution

The only things undetermined are C,k and y. The purpose is to find C for a given k and y. However I'm not too sure on the procedure in solving an eigenvalue expression like this so I tried to follow the Chebyshev method in Mathematica to solve for the eigenvalue C and I've ended up with a large output at the
Eigenvalues[result] part of the code.
that is too large to post full so here is a small part of it showing values and the general form:
$Root[...] \ , \ Root[-600. C^5 E^((3 y^2)/2000000000000) k^12 + <br /> 150000. C^4 E^((3 y^2)/2500000000000) k^12 + ...] \ ,Root[...]$

But I'm not sure whether I'm on the right track and if I am, how to go about plotting C vs k.

Thanks

Sorry if I'm not sounding clear

 

[primalfido]

Just thought I'd post the Mathematica code I'm working with:

Uo = 100
b = 6*10^(-11)
R = 60268000
r = R*Sin[76]
U = Uo*Exp[-b*(y^2)/(2*Uo)]
Ld = N[Sqrt[2*Uo/b]]
\[CapitalOmega] = N[(2*\[Pi]/(39 + 60*39 + 10*60*60))]
\[Beta] = 2*\[CapitalOmega]*Cos[(76) Degree]/R
Clear[y]

Uyy2 =  (b^{2} y^{2}/Uo)*Exp[-b*(y^2)/(2*Uo)] - b*Exp[-b*(y^2)/(2*Uo)]

DE = \[Phi]''[
   y] + (- k^2 +  (\[Beta] - D[U, {y, 2}] )/(U - g)) \[Phi][y]


t1 = DE /. {\[Phi][y] -> a[i]*ChebyshevT[i, y], 
   D[\[Phi][y], {y, a1_}] -> D[a[i]*ChebyshevT[i, y], {y, a1}]}

t2 = ExpandAll[t1]

n = 12

t3 = Sum[t2, {i, 0, n}]
t4 = Simplify[t3]
t5 = (t4 *(1/Sqrt[1 - y^2])*ChebyshevT[i, y] 2)/(\[Pi] c[i])
t6 = Table[t5, {i, 0, n, 2}];
t7 = (t6 /. c[0] -> 2) /. Table[c[i] -> 1, {i, 2, n}];
t8 = Collect[Expand[Numerator[t7]], y];
t9 = Transpose[Table[Coefficient[t8, y, i], {i, 0, 2 n, 2}]]
t10 = Table[y^i/(\[Pi] Sqrt[1 - y^2]), {i, 0, 2 n, 2}]
Timing[t11 = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \(1\)]\(t10 \
\[DifferentialD]y\)\)]
t12 = Simplify[Expand[t9.t11]]
t13 = Transpose[Table[Coefficient[t12, a[i]], {i, 0, n, 2}]]
b1 = Table[1, {i, 0, n, 2}]
t14 = Append[Table[t13[[i]], {i, 1, n/2}], b1]
s1 = Solve[t14[[n/2 + 1]].Table[a[i], {i, 0, n, 2}] == 0, a[n]]
t15 = Table[t14[[i]].Table[a[i], {i, 0, n, 2}], {i, 1, n/2}]
t16 = Expand[t15 /. s1]
t17 = Transpose[
  Table[Table[Coefficient[t16[[1, i]], a[j]], {i, 1, n/2}], {j, 0, 
    n - 2, 2}]]
(*How would I go about geting the complex phase speed g from this?*)
(*Like this? Am I approaching it wrong*)
vals = Eigenvalues[-t17 /. g -> 0]

I'm just not very sure on how to go about solving this complex eigenvalue problem numerically.