Very different results using filter object and b,a coefficients

I'm comparing the output of digital filtering using MATLAB filter object and b-a coefficients using tf function, and they are really different. I'll appreciate if anyone could point me out in the right direction. Here there are my results and the test script (download data.txt):

srate=64;
freqrange=[0.4 3.5];
file='data.txt';
load(file);
file=split(file,'.');
file=file{1};
data=eval(file);

st=srate*60;
ed=srate*60*2;
data=data(st:ed);%1 minute data
m=numel(data);

x=data;
R=0.1;%10% of signal
Nr=50;
NR=min(round(m*R),Nr);%At most 50 points
x1=2*x(1)-flipud(x(2:NR+1));%maintain continuity in level and slope
x2=2*x(end)-flipud(x(end-NR:end-1));
x=[x1;x;x2];

[xx,dbp]=bandpass(x,freqrange,srate);%same result if use filter(dbp,xx)
data_fil=xx(NR+1:end-NR);

[b,a]=dbp.tf;
xx=filter(b,a,x);
data_fil_ba=xx(NR+1:end-NR);

f=figure;
s(1)=subplot(2,1,1,'Parent',f);
plot(s(1),[data data_fil])
title(s(1),'dbp')
s(2)=subplot(2,1,2,'Parent',f);
plot(s(2),[data data_fil_ba])
title(s(2),'ba')
linkaxes(s,'x');

and here there is the result I obtain:

This is caused by numerical instability in the [b,a] filter representation, from the Matlab docs:

Numerical Instability of Transfer Function Syntax

In general, use the [z,p,k] syntax to design IIR filters. To analyze or implement your filter, you can then use the [z,p,k] output with zp2sos. If you design the filter using the [b,a] syntax, you might encounter numerical problems. These problems are due to round-off errors and can occur for n as low as 4.

Looking at your filter with fvtool, you will see that the first plot (of the digitalFilter object) is smooth, where as the second plot, generated with the transfer function representation is very unstable in the lower frequency range. The magnitude of this filter even gets above 0 dB, meaning that the filter is unstable and will amplify these frequencies. That is why your results blows up, it goes to ~10^92.

The solution is to simply not use the transfer function ([b,a]) representation, and instead use the digitalFilter, or [z,p,k] representation. See here for more info on why to be careful on switching filter/model representations.

Code

srate=64;
freqrange=[0.4 3.5];
[~,dbp]=bandpass(rand(100,1),freqrange,srate);
[b,a]=dbp.tf;
fvtool(dbp)
fvtool(b,a)