Actually if you need the frequency where the gain takes the maximum, then things get a bit hairy. If it's a peaking EQ (ie. poles and zeroes at the same frequency) or a BLT band-pass/notch, then you will find the maximum either at pole frequency (boost) or DC/Nyquist (cut). For any other placement of zeroes, the actual maximum is typically NOT at the pole frequency, but rather somewhere "nearby" (and away from the zeroes) unless the Q is low enough that the maximum is found at DC or Nyquist.
For example, any biquad BLT low-pass (ie. both zeroes at -1) with Q>sqrt(.5) will have a maximum gain point higher than the DC gain somewhere in it's passband, yet the cutoff gain is Q, so it only reaches unity at Q=1. As Q increases, the maximum gain point gets closer to cutoff, but it only becomes exactly the same once you reach the self-oscillation limit at Q=inf.
As for the best method to actually compute the exact answer... I'm not sure, since this is rarely needed for audio (where you usually want to normalise by some other condition). You could probably figure out an expression for the magnitude response and take derivative with respect to w, then solve for the roots, but I'm not sure if this is workable this is in practice without numerical root finding. You could also do an inverse BLT first (which just warps frequencies, which is fine since the answer can always be warped back) to avoid some of the trigonometry hell, but again I'm not sure how much it would help.