As much as I like elliptic filters, creating three of them and filtering them in parallel is simply not efficient. Multiband filters are much easier to implement as FIR filters.
Try a FIR filter such as this one:
Fs = 100; % Sampling Frequency (Hz)
fcuts = [3.5 4 7 7.5 9.5 10 13 13.5 29 30 35 36]; % Frequencies
mags = [0 1 0 1 0 1 0 ]; % Passbands & Stopbands
devs = [0.05 0.01 0.05 0.01 0.05 0.01 0.05]; % Tolerances
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs); % Kaiser Window FIR Specification
n = n + rem(n,2);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); % Filter Realisation
figure
freqz(hh,1,2^14,Fs)
set(subplot(2,1,1), 'XLim',[0 50]); % Zoom Frequency Axis
set(subplot(2,1,2), 'XLim',[0 50]); % Zoom Frequency Axis
Then use the filtfilt function to do the actual filtering:
signal_filtered = filtfilt(hh, 1, signal);
The Bode plot for this filter:
