Two layer internal waves

This analysis is similar to a single layer system,  except there are 2 more unknowns and 2 more boundary conditions.  Properities or variables in the lower layer are denoted with a "2" and in the upper layer with a "1".

phi1:=f1*phase

phi2:=f2*phase

eta:=C*phase

phase=Exp[I(k*x-ω*t)];

f1=A1*Exp[k*z]+B1*Exp[-k*z];

f2=A2*Exp[k*z]+B2*Exp[-k*z];

bottombc=(D[f2,z]==0)/.z->-H2

topbc=(D[f1,z]==0)/.z->H1

dynbc=(ρ1 D[phi1,t]-ρ2 D[phi2,t]==-g*eta(ρ1-ρ2))/.z->0

kinbc1=(D[phi1,z]==D[eta,t])/.z->0

kinbc2=(D[phi2,z]==D[eta,t])/.z->0

oneeq=Simplify[Eliminate[{bottombc,dynbc,kinbc1,kinbc2,topbc},{A1,B1,A2,B2}]]

dis1=Solve[oneeq,ω]

Let's beat this into a nice dispersion relation:

dis2=ω^2==(ω^2/.dis1[[1,1]])

dis3=FullSimplify[dis2]

dis4=dis3//ExpToTrig//FullSimplify

Now just for fun let's look at a limit that brings us back to a some limits that may or may not look familar to you. (They looked familiar to me!)

dis5=dis4/.H2->H1

Limit[Csch[2 x]Sinh[x]^2//TrigToExp,x->Infinity]

So the limit is (ω^2)==(( g k (ρ2-ρ1) )/(ρ1+ρ2)).  Take a further limit ρ1⟶0, and it should look familar to you too!