DSolve leads to an equation that is different from paper solution

I'm following the method in the two images in an attempt to obtain a value for $\beta$. I know that there's a mistake in (4.36) from here so I used Mathematica's solution to continue the solution. I used the following code:

ClearAll[Y, y, \[Mu] , k, a, \[Theta]]
sol = DSolve[{Y''[y] - ((\[Mu] k \[Pi])/a)^2 Y[y] == (-8 \[Theta])/(
     k \[Pi]), Y[-b/2] == 0, Y[b/2] == 0}, Y, y];
Y2[y_] = FullSimplify[ExpToTrig[Y[y] /. sol[[1]]]];
\[Phi][x_, y_] = Y2[y]*Sin[(k \[Pi])/a x];
gj = 2 \!\(
\*SubsuperscriptBox[\(\[Integral]\), 

FractionBox[\(-b\), \(2\)], 
FractionBox[\(b\), \(2\)]]\(\((
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(a\)]\[Phi][x, 
         y] \[DifferentialD]x)\) \[DifferentialD]y\)\);
a = b; Gx = Gy; \[Mu] = Sqrt[Gx/Gy];
\[Beta] = gj/(Gx a b^3)

to get that

$$\beta=\frac{32 \theta \sin ^2\left(\frac{\pi k}{2}\right) \left(\pi b k-2 b \tanh \left(\frac{\pi k}{2}\right)\right)}{\pi ^5 b G_y k^5}$$ According to the solution in the images and my understanding, $\beta$ is a a factor so setting $a=b$ and $G_y=G_x$ should lead to an elimination of those terms, i.e. they divide and become 1 like in the case of the $c$ factor in the paper, but that is not the case in the solution using Mathematica. There are an extra $b$, $\theta$, and $G_y$ that cannot be cancelled with $a$ and $G_x$.

When the author of the paper does N[Sum[\[Beta], {k, 1, 60, 2}]] (they didn't use Mathematica but it's to get the point across) setting $c=1$ they the following result

which agrees with simulations. When I run N[Sum[\[Beta], {k, 1, 60, 2}]] I get $\frac{0.140577}{G_y}$ which will give me a different answer to the author of the paper.

This is the 4th time that this happens. I've been using several techniques and following various similar solutions but I can't seem to get the equations from Mathematica and from the papers to agree. Is there something wrong with my code that I'm not seeing? Can I rearrange my equation for $\beta$ so that I won't end up with an extra $G_y$ in the solution? Or is it simply that the author made a mistake?

The governing equation for Prandtl’s stress function $\phi(x,y)$ is

$$\phi_{xx}+\frac{G_{zy}}{G_{zx}}\phi_{yy}=-2\theta$$ with BCs $\phi(\pm a/2,y)=\phi(x,\pm b/2)=0$.