Уважаемые эксперты! Пожалуйста, ответьте на вопрос:

Здравствуйте! У меня вопрос по поводу матлаба. Я не нашёл в разделах Матлаб. Делаю программу в Матлаб по матстатистике. Есть готовый алгоритм, частично реализован в коде, застрял на вычислениях с матрицами и итерациями с ними. Помогите, пожалуйста. Благодарю за внимание!

Приложение:

[fname,pname] = uigetfile('D:\MATL\QED60.txt','Открыть файл');
name = [pname fname];
b = load(name);
N = length(b);
dvs = length(b(1,:));
g = (b(1:N, 2:dvs))';
dvst = dvs - 1;
int = 1;
hint = int / dvst;
x=hint:hint:int;
[fname,pname] = uigetfile('D:\MATL\60 РД.txt','Открыть файл');
name = [pname fname];
b1 = dlmread(name);
dvs1 = length(b1(1,:));
g1 = (b1(1:N, 2:dvs1))';
dvst1 = dvs1 - 1;
hint1 = int / dvst1;
x1=hint1:hint1:int;
n = 7;m = 4;
k = 1;
p = figure;
grid
for i=1:N
clc;
y = g(:,i);
r = g1(:,i);

A = [1 0.142857143 0.0204081633 0.00291545;
1 0.285714286 0.0816326532 0.023323615;
1 0.428571429 0.1836734697 0.078717201;
1 0.571428572 0.3265306128 0.186588921;
1 0.714285715 0.5102040826 0.364431487;
1 0.857142858 0.7346938790 0.629737611;
1 1 1 1 ];
U1 = inv(A'*A) * A'*y; %также проблема в итерации U, вопрос мой в том, что
U11 = U1(1); %возможна ли итерация матриц?
U12 = U1(2);
U13 = U1(3);
U14 = U1(4);
ydet = (U11+U12*x+U13*x.^2+U14*x.^3)';
sig = sum((ydet-y).^2)/(n-m) ; %рассчёт первой дисперсии
if abs(y-ydet) <= k*sig^0.5
y11= y;
v(j,1)= 1
elseif y-ydet > k*sig^0.5
y11(j,1) = ydet(j) + (k*sig^0.5);
v(j,1) = 0
else y-ydet < -(k*sig^0.5)
y11 = ydet - (k*sig^0.5);
v(j,1) = 0
end
y1 = v(j,1) * y11'; %проблема здесь с вычислением матрицы. нужна матрица 7х1
koln = sum(v(:,1) == 1)

U2 = inv(AtA) * At*y1;
U21 = U2(1);
U22 = U2(2);
U23 = U2(3);
U24 = U2(4);
ydet2 = (U(21) + U(22)*x + U(23)*x.^2 + U(24)*x.^3)';

sig2 = sum((ydet2-y1).^2)/(koln-m);
for j=1:n
if abs(y1-ydet2)<= k*sig2^0.5
y2 = y;
elseif y1-ydet2 > k*sig2^0.5
y2 = ydet2 + (k*sig2^0.5);
else I21 < -(k*sig2^0.5)
y2 = ydet2 - (k*sig2^0.5);
end
end
U3 = inv(AtA) * At*y2;
U31 = U3(1);
U32 = U3(2);
U33 = U3(3);
U34 = U3(4);
ydet3 = (U31 + U32*x + U33*x.^2 + U34*x.^3)';

cla
hold on;
lin=plot(x,y,'-+',x1,r,'-*');
set(lin,'linewidth',1.5);
plot(x,oblp,x,oblm)
plot(x,ydet,x,ydet2,x,ydet3,x)
xlabel('Высота, м')
ylabel('Относительная мощность, МВт')
l=legend('Показания ДПЗ','Расчётные данные','Верхняя граница','Нижняя граница','ydet','ydet2','ydet3','Location','NorthEastOutside');
title(l,'Определения')
eval(['print -dpng fig' int2str(i)]); %сохраняет графики в файлах
end
hold off