Построить пилот, по заданным данным

Сроду не программировал и сразу кинули в matlab. Не знаю в чем здесь проблема. Буду рад если поможете. Извиняюсь за такое кол-во строк кода :(

clear all
clc
NumVar = 10;
%1. Зададим исходные данные:
%Геометрические параметры пластины:
Lx = 0.1*NumVar;
Ly = 0.2*NumVar;
Lz = 0.04*NumVar;
% Слои пластины:
z1 = 0.6*Lz;
z2 = 0.7*Lz;
z_zv = 0.8*Lz;
%Коэффициент температуропроводности:
a1 = NumVar*1*10^-5;
a2 = NumVar*2*10^-4;
a3 = 5*a1;
% Дискретизация:
Nx = 5;
Ny = 5;
Nz = 15;
delX = Lx/(Nx-1);
delY = Ly/(Ny-1);
delZ = Lz/(Nz-1);
delZ1 = (6/15)*(Lz/5);
delZ2 = (3/15)*(Lz/2);
delZ3 = (6/15)*(Lz/5);
delt = 0.01;
%Технические параметры:
lambtp1 = NumVar*0.02;
lambtp2 = NumVar*5;
t=0;
i=3;
Fold = zeros(Nx,Ny);
T_g = zeros(Nx,Ny);
lapl = zeros(Nx,Ny);
T = zeros(Nx,Ny,Nz);
delT_etta_gamma_ksi = zeros(Nx,Ny,Nz);
U = zeros(Nx,Ny);
%2. Расчет данных:
%Расчет управления:
N1 = 1;
N4 = 1509.977;
E1 = 1.39;
E2 = 934.1678;
E4 = 0.00099;
N1 = 80.765158;
N4 = 57.231289;
E1 = 34.8664826;
E2 = 9.278370;
E4 = 0.05944;
SUMM = zeros(5,5);
SUMM_prev = zeros(5,5);
for lmb = 1:1*10^5
t = t + delt;
for etta = 2:Nx-1
for gamma = 2:Ny-1
T_g(etta,gamma) = 10;
end
for etta = 1:Nx
for gamma = 1:Ny
F(etta,gamma) = T_g(etta,gamma) - T(etta,gamma,4)
end
for etta = 2:Nx-1
for gamma = 2:Ny-1
lapl(etta,gamma) = (F(etta-1,gamma)-2*F(etta,gamma)+F(etta+1,gamma))/delX^2+...
(F(etta,gamma-1)-2*F(etta,gamma)+F(etta,gamma+1))/delY^2;
SUMM_prev(etta,gamma) = E4*(((N4-1)/N4)*F(etta,gamma)-(1/N4)*lapl(etta,gamma))*delt;
SUMM(etta,gamma) = SUMM(etta,gamma)+SUMM_prev(etta,gamma);
% U(etta,gamma) = 1000*cos(pi*(i-0.5)*etta/(Nx-1))*sin(pi*gamma*(i-0.5)/(Ny-1));
U(etta,gamma) =E1*(((N1-1)/N1)*F(etta,gamma)-(1/N1)*lapl(etta,gamma)) + SUMM(etta,gamma) + E2*((F(etta,gamma)-Fold(etta,gamma))/delt);
end
%Тепловые процессы, протекающие в печи
for etta = 2:Nx-1
for gamma = 2:Ny-1
for ksi = 2:6
delT_etta_gamma_ksi(etta,gamma,ksi) = a1*delt*((T(etta-1,gamma,ksi)-2*T(etta,gamma,ksi)+T(etta+1,gamma,ksi))/delX^2 + ...
(T(etta,gamma-1,ksi)-2*T(etta,gamma,ksi)+T(etta,gamma+1,ksi))/delY^2 + .
(T(etta,gamma,ksi-1)-2*T(etta,gamma,ksi)+T(etta,gamma,ksi+1))/delZ1^2);
end
for etta = 2:Nx-1
for gamma = 2:Ny-1
for ksi = 8
delT_etta_gamma_ksi(etta,gamma,ksi) = a2*delt*((T(etta-1,gamma,ksi)-2*T(etta,gamma,ksi)+T(etta+1,gamma,ksi))/delX^2 + ...
(T(etta,gamma-1,ksi)-2*T(etta,gamma,ksi)+T(etta,gamma+1,ksi))/delY^2 +
(T(etta,gamma,ksi-1)-2*T(etta,gamma,ksi)+T(etta,gamma,ksi+1))/delZ2^2);
end
for etta = 2:Nx-1
for gamma = 2:Ny-1
for ksi = 10:14
delT_etta_gamma_ksi(etta,gamma,ksi) = a3*delt*((T(etta-1,gamma,ksi)-2*T(etta,gamma,ksi)+T(etta+1,gamma,ksi))/delX^2 + ...
(T(etta,gamma-1,ksi)-2*T(etta,gamma,ksi)+T(etta,gamma+1,ksi))/delY^2 +
(T(etta,gamma,ksi-1)-2*T(etta,gamma,ksi)+T(etta,gamma,ksi+1))/delZ3^2);
end
%Граничные условия при переходе из среды в среду
for etta = 2:Nx-1
for gamma = 2:Ny-1
T(etta,gamma,7) = (lambtp1*T(etta,gamma,6)+lambtp2*T(etta,gamma,8))/(lambtp1+lambtp2);
T(etta,gamma,9) = (lambtp1*T(etta,gamma,10)+lambtp2*T(etta,gamma,8))/(lambtp1+lambtp2);
end
for ksi = 2:Nz-1
for gamma = 2:Ny-1
for etta = 2:Nx-1
T(etta,gamma,ksi) = T(etta,gamma,ksi)+delT_etta_gamma_ksi(etta,gamma,k si);
end
%Входные воздействия
for etta = 2:Nx-1
for gamma = 2:Ny-1
% T(etta,gamma,1) = delZ*utp/lambtp1 + T(etta,gamma,2);
T(etta,gamma,1) = delZ*U(etta,gamma)/lambtp1 + T(etta,gamma,2);
end
%Граничные условия:
for etta = 1:Nx
for gamma =1:Ny
for ksi = 1:Nz
T(Nx,gamma,ksi) = T(Nx-1,gamma,ksi);
T(etta,Ny,ksi) = T(etta,Ny-1,ksi);
T(etta,gamma,Nz) = T(etta,gamma,Nz-1);
F(Nx,gamma) = T(Nx-1,gamma);
F(etta,Ny) = T(etta,Ny-1);
end
for etta = 1:Nx
for gamma = 1:Ny
Fold(etta,gamma) = F(etta,gamma);
end
end
ForPlotT(lmb, = T(4,2,6);
figure(1), plot(lmb,T(4,4,12)), hold on
figure(2), plot(lmb,U), hold on
end
%Расчет характеристик:
n = length(ForPlotT);
M = ForPlotT(n);
M_vv = U(4,4);
K = M/M_vv
Psi_x = pi*(i-0.5)/Lx;
Psi_y = pi*(i-0.5)/Ly;
G = Psi_x^2 + Psi_y^2