Matlab - Guardar matrices dentro de otras de forma iterativa

casos=10;

liminf_p=0.3;

liminf_q=0.3;

limsup_p=0.9;

limsup_q=0.9;

for k=1:casos

P = liminf_p:((limsup_p-liminf_p)/casos):limsup_p;

Q = liminf_q:((limsup_q-liminf_q)/casos):limsup_q;

Sinf(k)= P(k)+Q(k)*i;

IL(k) =Sinf(k)/Vinf;

Vg(k)= Vinf+IL(k)*Ze;

Vt_ang_rad(k)= angle(Vg(k));

Vt_mod(k)= abs(Vg(k));

Ig(k)=(cos(Vt_ang_rad(k))+i*sin(Vt_ang_rad(k))-Vbinf_mod)/(Re+i*Xe);

E(k)=cos(Vt_ang_rad(k))+i*sin(Vt_ang_rad(k))+Ig(k)*(Rs+i*Xq);

ang_pot(k)=angle(E(k));

ang_pot_rad(k) = (ang_pot(k)*2*pi/360);

%tensiones y corrientes en ejes dq

Id_Iq(k)=(cos((pi/2)-ang_pot(k))+i*sin((pi/2)-ang_pot(k)))*Ig(k);

Vd_Vq(k) = (cos((pi/2)-ang_pot(k))+i*sin((pi/2)-ang_pot(k)))*(cos(Vt_ang_rad(k))+i*sin(Vt_ang_rad(k)));

% Eq transiente

Eq_trans(k) = imag(Vd_Vq(k))+Xd_trans*real(Id_Iq(k));

%Efd

Efd(k)=Eq_trans(k)+(Xd-Xd_trans)*real(Id_Iq(k));

% Vref

Vref(k) = Vt_mod(k)+(Efd(k)/Ka);

%Tm

Tm(k)= imag(Id_Iq(k))*(Eq_trans(k)+(Xq-Xd_trans)*real(Id_Iq(k)));

delta= sqrt(Re)+(Xe+Xq)*(Xe+Xd_trans);

% Constantes

K1(k) = -(imag(Id_Iq(k))*Vbinf_mod*(Xd_trans-Xq)*((Xq+Xe)*sin(ang_pot(k))-Re*cos(ang_pot(k)))+Vbinf_mod*((Xd_trans-Xq)*real(Id_Iq(k))-Eq_trans(k))*((Xd_trans+Xe)*cos(ang_pot(k))-Re*sin(ang_pot(k))))/delta;

K2(k) = (imag(Id_Iq(k))*delta-imag(Id_Iq(k))*(Xd_trans-Xq)*(Xq+Xe)-Re*(Xd_trans-Xq)*real(Id_Iq(k))+Re*Eq_trans(k))/delta;

K3_inv(k)=1+((Xd-Xd_trans)*(Xq+Xe))/delta;

K4(k) =((Xq+Xe)*sin(ang_pot(k))-Re*cos(ang_pot(k)))*Vbinf_mod*(Xd-Xd_trans)/delta;

K5(k) = ((real(Vd_Vq(k))/Vt_mod(k))*Xq*(Re*Vbinf_mod*sin(ang_pot(k))+Vbinf_mod*cos(ang_pot(k))*(Xd_trans+Xe))+(imag(Vd_Vq(k))/Vt_mod(k))*Xd_trans*( Re*Vbinf_mod*cos(ang_pot(k))-Vbinf_mod*sin(ang_pot(k))*(Xe+Xq)))/delta;

K6(k)=(((real(Vd_Vq(k))/Vt_mod(k))*Xq*Re-(imag(Vd_Vq(k))/Vt_mod(k))*Xd_trans*(Xe+Xq))/delta)+(imag(Vd_Vq(k))/Vt_mod(k));

K3(k)= (K3_inv(k))^(-1);

% matriz barra infinita AVR sin considerar PSS

A_patron = [(-d*Ws/(2*H)) (-K1(k)/(2*H)) (-K2(k)/(2*H)) 0;Ws 0 0 0;0 (-K4(k)/Tdo_trans) (-K3_inv(k))*(1/Tdo_trans) (1/Tdo_trans);0 ((-Ka*K5(k))/Ta) ((-Ka*K6(k))/Ta) (-1/Ta)];

%matrices para crear las tf para graficar los eig de los casos

B=[0;0;0;Ka/Ta];

C=[1 0 0 0];

D=zeros(1);

sys=[];

ay=0.1;

by=7.2079;

gy=26.3157;

Kpss= 0.000104;

Ac=[-ay 0 0 0;gy-ay*by -gy 0 0;by*gy-ay*by^2 gy-by*gy -gy 0;by^2*gy-ay*by^2 by*gy-by^2*gy gy-by*gy -gy];

Bc=[1 by by^2 by^3]';

Cc=[0 0 0 Kpss];

Ax=[];

for u=1:casos

Ax=[Ax A_patron]; % matriz que guarda la matriz A_patron para cada caso

end

% EN ESTA PARTE PRESENTO EL PROBLEMA MENCIONADO

Am=[];

for j =1:casos 

Am(j,j)=[Ax(:,4*j-3:4*j) B*Cc;Bc*C*Ax(:,4*j-3:4*j)  Ac];

end

end