Ускорить вычисление матрицы

Я работаю над линейной моделью прогностического управления, и мне нужно вычислить некоторые матрицы только для контроллера.. вычисление одного из них занимает много времени, и я хотел бы спросить, есть ли лучший способ кодирования этого вычисления. Я использую MATLAB, но я также понимаю FORTRAN.

Ну, я хочу вычислить матрицу (Φ), но способ, которым я это делаю, требует много времени, чтобы вычислить ее. Матрица Φ имеет вид (правильный):
MPC_matrices.
здесь есть книга, где я нашел это изображение в случае, если вам нужно обратиться (особенно Страница 8).

Теперь код, который я написал в MATLAB, выглядит так: (переместил его после EDIT)

Учитывая, что у меня будут довольно большие переменные NS, Np и Nc, это займет много времени, чтобы сделать этот расчет. Есть ли оптимальный способ (или, по крайней мере, лучший, чем мой) ускорить этот расчет?

EDIT

После рассмотрения комментариев @Daniel's & user2682877 я протестировал это

clear;clc

Np = 80;

Nc = Np / 2;

m  = 3;

q  = 1;

Niter = 30;

MAT = zeros(Niter,5);

for I=1:Niter

    NS = 10 * I;

    A = rand(NS,NS);

    B = rand(NS,m);

    C = rand(1,NS);

    tic

    Phi1 = zeros(Np*q,Nc*m);

    CB = C * B;

    for i=1:Np

        for j=1:Nc

            if j<i

                Phi1( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = C * A^(i-1-(j-1)) * B;

            elseif j==i

                Phi1( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = CB;

            end

        end

    end

    t1 = toc;

% предложение Даниэля

    tic

    Phi2=zeros(Np*q,Nc*m);

    CB = C * B;

    for diffij=0:Np-1

        if diffij>0

            F=C * A^diffij * B;

        else

            F=CB;

        end

        for i=max(1,diffij+1):min(Np,Nc+diffij)

            j=i-diffij;

            Phi2( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = F;

        end

    end

    t2 = toc;

% user2682877 предложение

    tic

    Phi3=zeros(Np*q,Nc*m);

    temp = B;

    % 1st column

    Phi3( (q*1-(q-1)):(q*1) , (m*1-(m-1)):(m*1) ) = C * B;

    for i=2:Np

        % reusing temp

        temp = A * temp;

        Phi3( (q*i-(q-1)):(q*i) , (m*1-(m-1)):(m*1) ) = C * temp;

    end

    % remaining columns

    for j=2:Nc

        for i=j:Nc

            Phi3( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) =...

                Phi3( (q*(i-j+1)-(q-1)):(q*(i-j+1)) , (m*1-(m-1)):(m*1) );

        end

    end

    t3 = toc;



    MAT(I,:) = [I, NS, t1, t2 ,t3];

    fprintf('I iteration = %g n', I);

end

figure(1)

clf(1)

hold on

plot(MAT(:,2),MAT(:,3),'b')

plot(MAT(:,2),MAT(:,4),'r')

plot(MAT(:,2),MAT(:,5),'g')

hold off

legend('My <Unfortunate> Idea','Daniel`s suggestion','user2682877 suggestion')

xlabel('NS variable')

ylabel('Time, s')

И вот результат::

Имейте в виду, что теперь NS = 300 ,но по мере того, как я расширяю свою модель (я намерен включать все больше и больше уравнений и переменных в модель пространства состояний), эти переменные (в основном NS и Np) будут все больше и больше.

@ Daniel's 2nd comment, я знаю, что выполняю больше вычислений, чем следовало бы, но у меня нет опыта ограничивает мои представления о том, как это сделать.

комментарий @ durasm, я не совсем знаком с parfor, но я проверю его.

Ссылаясь на ответы: я проверю ваши предложения, как только я их пойму (...) и вернуться к вам.

Результаты
Очевидно, что моя первоначальная мысль лишь немного хуже той, что была высказана здесь.. Спасибо вам, ребята! Вы были очень полезны!