Electronic – Least Square Fit – Transfer Function

I have here this continuous transfer function:

$$\frac{Y(s)}{U(s)} = \frac{s+3}{s^3 + 12 s^2 + 39 s + 28}$$

Passing this tf to discrete domaing like this:

T = 0.01;
sys_discrete = c2d(sys, T);

I got this:

$$\frac{Y(z)}{U(z)} = \frac{4.853.10^{-5} z^2 + 1.57.10^{-11} z – 4.57.10^{-5}}{z^3 – 2.883 z^2 + 2.77 z – 0.8869}$$

Considering this:

$$\frac{Y(z)}{U(z)} = \frac{a_0 z^2 + a_1 z + a_2}{b_0 z^3 +b_1 z^2 + b_2 z + b_3}$$

and doing the inverse Z transform, I got this:

$$y(k) = \frac{a_0}{b_3}u(k-2) + \frac{a_1}{b_3}u(k-1) + \frac{a_2}{b_3}u(k) – \frac{b_0}{b_3}y(k-3) – \frac{b_1}{b_3}y(k-2) – \frac{b_2}{b_3}y(k-1) \\
y(k) = \begin{bmatrix}u(k-2) & u(k-1) & u(k) & y(k-3) & y(k-2) & y(k-1)\end{bmatrix}.\begin{bmatrix} a_0/b_3 \\ a_1/b_3 \\ a_2/b_3 \\ -b_0/b_3 \\ -b_1/b_3 \\ -b_2/b_3 \end{bmatrix}
$$

Renaming these guys:

$$ \begin{bmatrix} a_0/b_3 \\ a_1/b_3 \\ a_2/b_3 \\ -b_0/b_3 \\ -b_1/b_3 \\ -b_2/b_3 \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \\ c_5 \\ c_6 \end{bmatrix} $$

So when I find this C vector, I will back in this way (just replacing back):

$$\frac{Y(z)}{U(z)} = \frac{c_1 z^2 + c_2 z +c_3}{-c_4 z^3 -c_5 z^2 -c_6 z + 1}$$

So what I'm trying is to find C:

$$Y = A.C \\ A^T.Y = A^T.A.C \\ (A^T.A)^{-1}.(A^T.Y) = (A^T.A)^{-1}.(A^T.A).C \\ (A^T.A)^{-1}.(A^T.Y) = I.C \\ (A^T.A)^{-1}.(A^T.Y) = C$$

That's how I'm doing this at MATLAB:

t = 0:T:10; %time vector
u = wgn(length(t), 1, 1); %input vector
[y, t] = lsim(sys, u, t); %output vector applied to sys

%[u(k-2) u(k-1) u(k) y(k-3) y(k-2) y(k-1)]
A = [u(2:end-2) u(3:end-1) u(4:end) y(1:end-3) y(2:end-2) y(3:end-1)]; %doing the right shifts here, no problem about losing data...
Y = y(4:end);
C=(A'*A)\(A'*Y); %finding C
c1 = C(1); c2 = C(2); c3 = C(3);
c4 = C(4); c5 = C(5); c6 = C(6);
sys_ls = tf([c1 c2 c3], [-c4 -c5 -c6 1], T);

When I step(sys, sys_ls), that's what I'm getting:

I was hoping to find something very very similar… What am I missing to get this system fit right?