Using IIR filters to do that requires connecting all three in series, and a low-order Butterworth design is likely not going to meet your expectations, since the frequencies are ao close together.
I would use a single FIR filter instead:
Fs = 250;
fcomb = [[49 49.5 50.5 51], [49 49.5 50.5 51]+10, [49 49.5 50.5 51]+12.5];
mags = [[1 0 1], [0 1], [0 1]];
dev = [[0.5 0.1 0.5], [0.1 0.5], [0.1 0.5]];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh, 1, 2^20, Fs)
set(subplot(2,1,1), 'XLim',[45 65])
set(subplot(2,1,2), 'XLim',[45 65])
Experiment with this to get different results. This designs a 420-order FIR filter, and decreasing the stopband frequency widths will increase its length. (Increasing them would cause them to overlap, and that will throw an error.) It will still be an efficient filtter, however, compared with dascading three IIR filters in a loop.
Use filtfilt to do the actual filtering, Use the same calling convention for the first two arguments as I use in the freqz call.