Идентификация объекта управления по методу наименьших квадрантов (2-й вариант)

%задаюсь значениями для трансфер функции объекта s1 ( это то что пишется в комадной строке матлаба )
num=25;
den=[ 35 15 1];
%задаюсь трансфер функции желаемой передаточной функции
wzh=tf([2],[4 2 2]);

[s2,wr] = mnk(num,den,wzh)

function [s2,wr] = mnk(num,den,wzh)
s1=tf(num,den);
T_end=60; % интервал измерений
dt=0.2; % шаг дискретизации
t=0:dt:T_end; % массив дискретного времени
N=length(t); % размер выборки
u=ones(N,1); % массив значений единичного входного воздействия
y=lsim(s1,u,t); %массив значений выходного воздействия
n = max([length(num)-1,length(den)-1]);
I=diag([1 1 1]);
i=1; % начальный шаг
P=1000*I; % начальное приближение
betta=[0;0;0];
bet(i,:)=betta; % массив оценок параметров
% формирование расширенной матрицы
for i=n:N-1
R=[y(i+n-2:-1:i-1);u(i+n-2:-1:i)]'; 
gamma=P*R'/(R*P*R'+1);
betta=betta+gamma*(y(i+1)-R*betta);
P=(I-gamma*R)*P;
bet(i,:)=betta;
end;
% расчет параметров непрерывной модели
T1=dt^2/(1-betta(1)-betta(2)); 
T2=(betta(2)*T1+T1+dt^2)/dt;
K=betta(3)*T1/dt^2;
s2=tf([K],[T1 T2 1]);%передаточная функция идентифицированного объекта
y2=lsim(s2,u,t);
wr=[wzh/s2];%передаточная функция регулятора
% сравнение переходных характеристик объекта и модели
plot(t,y,t,y2,':'); 
grid;

num=25;den=[36 1];

??? Subscript indices must either be real positive integers or
logicals.
 
Error in ==> mnk at 17
R=[y(i+n-2:-1:i-1);u(i+n-2:-1:i)]';

a=20