Расчет элементов матрицы и циклическая нумерация

for I=1:54
t=0;
U(I,1)= 0;
U(I,2)= 0;
for J = 2:10
if rem(I,3)~=0
q=sin(7*I*t);
else
q = 0;
U(I,J)=inv(M(I,I)) * (-M(I,I) * q - K(I,I) * U(I,J-1)) * t*t + E(I) * (2*U(I, J-1) - U(I, J-2));
t= t + 0.001;
end;
end;

clear, clc
 
M = rand(54);
K = rand(54);
E = @(t) 1;
for I=1:54
    t=0;
    U(I,1)= 0;
    U(I,2)= 0;
    for J = 3:10
        if rem(I,3)~=0
            q=sin(7*I*t);
        else
            q = 0;
        end;
        U(I,J)=inv(M(I,I)) * (-M(I,I) * q - K(I,I) * U(I,J-1)) * t*t + E(I) * (2*U(I, J-1) - U(I, J-2));
        t= t + 0.001;    
    end;
end
surf(U)

M = rand(54,54);
K = rand(54,54);
E = eye (54,54);
for I=1:54
    t=0;
    U(I,1)= 0;
    U(I,2)= 0;
    for J = 3:10000
        if rem(I,3)~=0
            q=sin(7*I*t);
        else
            q = 0;
        end;
        U(I,J)=inv(M(I,I)) * (-M(I,I) * q - K(I,I) * U(I,J-1)) * t*t + E(I,I) * (2*U(I, J-1) - U(I, J-2));
        t= t + 0.001;    
    end;
end

clear, clc
 
M = 5+rand(54);
K = diag(2*ones(1,54))+diag(ones(1,51),3)+diag(ones(1,51),-3);
E = eye(54);
 
n = 1e4; 
U = zeros(54,n); % заготовка из нулей
t = 0:0.001:(n-1)*0.001; % массив времени
for i = 2:n-1 
        q = repmat([sin(7*i*t(i));0;0],18,1); % повторяем вектор 18 раз
        U(:,i+1) = inv(M)*(-M*q - K*U(:,i))*t(i)^2 + E*(2*U(:,i) - U(:,i-1));
end