How to plot interaction diagram

Hello everyone,

I want to plot an interaction diagram for biaxially loaded column. It should be like in picture but shaded parts are not necessary. I have a script file to calculate N, Mx and My values and there is another question here that I want to record these values for every loop. There are 2 loops, first one take kyh and and calculate N,Mx,My in second loop t and its go on like that. So, for every different kyh and t set, I obtain different N,Mx,My set but I don't know how to record to a matrix or an excel file in these loops and how to plot shape like that.

Thanks for your help.

%Interaction Diagram of Biaxially Loaded Columns
clear
clc
%Inputs
h=300; %height of column section
b=300; %width of column section
fck=20; %characteristic strength of concrete
fyk=420; %characteristic strength of steel
ds=20; %diameter of steel bar
k1=0.85; %k1 value w.r.t concrete strength

%Location of bars
XY=zeros(4,2);
o=2;
for i=1:4
    XY(i,1)=T(o,1);
    XY(i,2)=T(o,2);
    o=o+1;
end

%Initial Settings
v=4*h;
p=zeros(4,1); %distances of steel bars from inclined n.a.
EpsS=zeros(4,1); %Strain of steel bars
SgmS=zeros(4,1); %Stress of steel bars
Fs=zeros(4,1); %Steel forces
Msx=zeros(4,1); %Steel moment about x
Msy=zeros(4,1); %Steel moment about y
EpsSy=fyk/200000;
As=(pi*ds^2)/4;
%Algorithm
for kyh=1:v
    for t=1:89
        c=kyh*(cosd(t));
        
        %Steel Moment Calculations
        for i=1:4
            x=XY(i,1);
            y=XY(i,2);
            p(i)=(-x*sind(t))+(y*cosd(t))-(b*sind(t)/2)+((kyh-h/2)*cosd(t));
            EpsS(i)=0.003*p(i)/c;
            if (abs(EpsS(i)))<=EpsSy
                SgmS(i)=EpsS(i)*200000;
            else if (EpsSy<abs(EpsS(i)))&&(abs(EpsS(i))<=0.1)
                    SgmS(i)=fyk*EpsS(i)/(abs(EpsS(i)));
                else
                    SgmS(i)=0;
                end
            end
            SgmS(i)=round(SgmS(i));
            Fs(i)=SgmS(i)*As/1e3;
            Fs(i)=round(Fs(i));
            Msx(i)=Fs(i)*y/1e3;
            Msy(i)=Fs(i)*x/1e3;
        end
        TFs=sum(Fs);
        TMsx=sum(Msx);
        TMsy=sum(Msy);
        
        
        %Concrete stress block calculations
        a=k1*kyh;
        if a<=h
            f=a*cotd(t);
            if f<=b
                A=a*f/2;
                cgx=-(b/2-f/3);
                cgy=h/2-a/3;
            else
                z=a-b*tand(t);
                A=b*(a+z)/2;
                cgx=-(b/2-b*((2*z+a)/(z+a))/3);
                cgy=h/2-(z^2+z*a+a^2)/(3*(z+a));
            end
        else
            m=a*cotd(t);
            if m<=b
                n=(a-h)*cotd(t);
                A=h*(m+n)/2;
                cgx=-(b/2-(m^2+m*n+n^2)/(3*(m+n)));
                cgy=h/2-h*((2*n+m)/(n+m))/3;
            else
                g=b*tand(t)-(a-h);
                l=b-(a-h)*cotd(t);
                A=b*h-l*g/2;
                cgx=-(((b^2)*h/2-(l^2)*g/6)/(b*h-l*g/2)-b/2);
                cgy=((h^2)*b/2-(g^2)*l/6)/(b*h-l*g/2)-h/2;
            end
        end
        Fc=0.85*fck*A/1e3;
        N=Fc+TFs;
        Mx=Fc*cgy+TMsx;
        My=Fc*cgx+TMsy;
    end
end