Second order section digital Bandpass
So it seems that the formulas from the textbook for the bandpass are not accurate. 

I used following approach to calculate second order sections given an analogue filter lowpass. As an example following code displays the coefficients for a 4th order bandpass, calculated from a second order lowpass protopye.


Paramters

fs = 1000;
fs1 = 200;
fs2 = 300;
w1 = 2*pi*fs1/fs;
w2 = 2*pi*fs2/fs;
Omega1 = 2*tan(w1/2);
Omega2 = 2*tan(w2/2);
syms s z Om2 Om1 a1 b1 A0
% apply LP to BP transformation
G_s_BP = subs(H_s, s, (s^2 + Om1*Om2)/(s*(Om2 - Om1)));

Start with analogue transfer function H_s

H_s = A0 / (1 + a1*s + b1*s^2);
% apply LP to BP transformation
G_s_BP = subs(H_s, s, (s^2 + Om1*Om2)/(s*(Om2 - Om1)));

Digital conversion with bilinear transformation

G_z_BP = subs(G_s_BP, s, 2*(( z - 1)/(z + 1)));
% order polynomial such that coefficients can be read out
coeffs_z = collect(G_z_BP, z)

Those coeffs_z are very similar to the ones given by the iowahills software
The band-pass filter coefficients then can be calculated with any analogue low-pass prototype, just enter numerical values on read out coefficients from coeffs_z