A = [0 0.3; 0.4 0.5];
B = [2 ; 1];
C = [-1 1];
D = 0;
sys = ss(A,B,C,D,1);
N = 10;
u = idinput(N+1,'rbs');
[y,t,x] = lsim(sys,u);

lhs = [x(2:N+1,:)' ; y(1:end-1,:)'];
rhs = [x(1:N,:)' ; u(1:N,:)']

%ABCD = lhs*rhs'*inv(rhs*rhs') % normal equation

ABCD = lhs / rhs  % better numerics