Digital filtering in matlab not giving expected results

I'm trying to design and test a simple digital high pass filter in matlab. I have two scripts: the first one is to design the filter, the second one is the implementation fo the recursion algorithm

1st script: desgn of the filter:

s=sym('s');
z=sym('z');
w=sym('w');
p=sym('p');
f=sym('f');
x=sym('x');
Pi=sym('Pi');
Ts=sym('Ts');
vc=sym('vc');

n=2;
Fspb=1/(s^2+sqrt(2)*s+1);   % simple low pass butterworth
Fs=subs(Fspb, s, 2*Pi*vc/s); % tranform to high pass butterworth with vc cutting freq
Fz=subs(Fs, s, 2/Ts*(z-1)/(z+1)); %bilinear transformation
Fp=subs(Fz, z, exp(Ts*p))
Fw=subs(Fp, p, 1i*w);
Ff=subs(Fw, w, 2*Pi*f);

vmax=5000;
vc_val=1000; %1000 Hz
vccont=1/Pi/Ts*tan(Pi*Ts*vc_val);
Ts_val=0.0001
fval=0:0.1:vmax;

pretty(expand(Fz))

[num, den]=numden(Fz);
cn=coeffs(num,z)
cd=coeffs(den,z)

it gives me the coefficients (for z) and the frequency response of the filter:

cn =     
[ 1, -2, 1]

cd =
[ Pi^2*Ts^2*vc^2 - 2^(1/2)*Pi*Ts*vc + 1, 2*Pi^2*Ts^2*vc^2 - 2, Pi^2*Ts^2*vc^2 + 2^(1/2)*Pi*Ts*vc + 1] 

here is the mfilter displayed with the freqz function:

and my second script, implementing the filter and testing on a simple sinus function:

Pi=3.14;
v=500;
Ts=0.0001;
x=0:Ts:100000*Ts;

y=sin(x*2*Pi*v); % sinus, freq=v
figure;
plot(x,y);

%% filter def
vc_val=1; %1000 Hz
vccont=1/Pi/Ts*tan(Pi*Ts*vc_val);

a=Pi^2*Ts^2*vccont^2;
b=sqrt(2)*Pi*Ts*vccont;

%% filtering
yf=zeros(1,size(y,2));

y1=0; y2=0; y3=0; y4=0; x1=0; x2=0; x3=0; x4=0;

for i=3:size(y,2)
     if (i>1)
         y1=yf(i-1);
         x1=xyi-1);
         if (i>2)
             y2=yf(i-2);
             x2=y(i-2);
         end
     end
    yf(i)=1/(a-b+1)*(-(2*a-2)*y1-(a+b+1)*y2+y(i)-2*x1+x2);
end

figure;
plot(x,yf);

But it doesn't give me the result i was expecting (left: the original sinus, right: the filtered result at 1/2 of the cutting frequency):

I checked that my filter is stable, I can't see what the problem may be. Have you got any hint?

thanks!

I found the problem, which is that the numden() function of matlab doesn't give the right numerators and denominators. thanks to the people who took time to help me!

the following code gives a wrong result:

z=sym('z');
Pi=sym('Pi');
Ts=sym('Ts');
vc=sym('vc');

Fz=1/((Pi^2*Ts^2*vc^2)/(z^2 - 2*z + 1) + (2^(1/2)*Pi*Ts*vc)/(z - 1) + (2*Pi^2*Ts^2*vc^2*z)/(z^2 - 2*z + 1) + (Pi^2*Ts^2*vc^2*z^2)/(z^2 - 2*z + 1) + (2^(1/2)*Pi*Ts*vc*z)/(z - 1) + 1)

[num, den]=numden(Fz);
cn=coeffs(num,z)
cd=coeffs(den,z)

matlab gives the coefficients:

cn =

[ 1, -2, 1]

cd =

[ Pi^2*Ts^2*vc^2 - 2^(1/2)*Pi*Ts*vc + 1, 2*Pi^2*Ts^2*vc^2 - 2, Pi^2*Ts^2*vc^2 + 2^(1/2)*Pi*Ts*vc + 1]

which are obvisoulsy wrong

thanks!