plotting is giving me multiple lines

Greetings,

I made some modifications to my code and I am trying to plot this, but it is giving me different lines. My intent is to plot the first row of "stress1f(:,k)" for the x axis and y axis z_1=[z(k) z(k+1)]; I tried to do it inside and outside the for loop, but matlab seems to graph the same points multiple times when i do it inside the for loop and outside it seems the same. any suggestions would be greatly appreciated it

My graph should look as the three graphs on the right when I graph the first, second, and third row against z_1=[z(k) z(k+1)];

close all
clear all
clc
%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio

%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31

%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter

%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));

%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));

% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
    Q12,Q22,0;
    0,0,Q66];


%%Laminate rotations in degrees
theta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];

%%Z-Values
z=(-length(theta)/2+(0:length(theta)))'*h;

% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:length(theta)
    m=cosd(theta(i));
    n=sind(theta(i));
    T=[m.^2,n.^2,2.*m.*n;
        n.^2,m.^2,-2.*m.*n;
        -m.*n,m.*n,m.^2-n.^2];
      
    % Global transformed reduced stiffness coefficients Q_bar
    Q_bar{i}=inv(T).*Q.*T;
end

Aij=0;
Bij=0;
Dij=0;

for j=1:length(theta)
    A_j=Q_bar{j}*(z(j+1)-z(j));
    Aj{j}=A_j;
    Aij=Aij+Aj{j};
    B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
    Bj{j}=B_j;
    Bij=Bij+Bj{j};
    D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
    Dj{j}=D_j;
    Dij=Dij+Dj{j};
end

A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;

%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];

%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];

%%Shear
e_x=sym('Epsilonx');
e_y=sym('Epsilony');
gamma_xy=sym('Gammaxy');
shear=[e_x;e_y; gamma_xy];

%%Bending Twist
k_x=sym('Kx');
k_y=sym('Ky');
k_xy=sym('Kxy');
bending=[k_x;k_y;k_xy];

%%Shear Extension Coupling
SEC=N==A*shear;
SEC_F=solve(SEC);
e_xf=vpa(SEC_F.Epsilonx);
e_yf=vpa(SEC_F.Epsilony);
gamma_xyf=(SEC_F.Gammaxy);
shearf=[e_xf;e_yf; gamma_xyf];

%%Bending
BEC=M==D*bending;
BEC_F=solve(BEC);
k_xf=vpa(BEC_F.Kx);
k_yf=vpa(BEC_F.Ky);
k_xyf=vpa(BEC_F.Kxy);
bendingf=[k_xf;k_yf;k_xyf];

for k=1:length(theta)
 stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf;%%BOTTOM LAYER
 stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf;%%TOP LAYER
 stress1f(:,k)=stress1; %Bottom
 stress2f(:,k)=stress2; %Top
 b=stress1f(1,:)';
 t=stress2f(1,:)';
 tt=[b t];
 z_1=[z(k) z(k+1)];
%  figure(1)
%  stressx=[stress1f(1,k) stress2f(1,k)];
%  z_1=[z(k) z(k+1)];
%  plot(stressx,z_1)
%  xline(0)
%  hold on 
%  figure(2)
%  stressy=[stress1f(2,k) stress2f(2,k)];
%  z_2=[z(k) z(k+1)];
%  plot(stressy,z_2)
%  xline(0)
%  figure(3)
%  stressxy=[stress1f(3,k) stress2f(3,k)];
%  z_3=[z(k) z(k+1)];
%  plot(stressxy,z_3)
%  xline(0)
%  hold on
end
figure(1)
plot(tt,z_1)

Answer 1

[Alan Stevens]

A little more like this perhaps:

%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio

%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31

%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter

%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));

%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));

% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
    Q12,Q22,0;
    0,0,Q66];


%%Laminate rotations in degrees
theta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];

%%Z-Values
z=(-length(theta)/2+(0:length(theta)))'*h;

% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:length(theta)
    m=cosd(theta(i));
    n=sind(theta(i));
    T=[m.^2,n.^2,2.*m.*n;
        n.^2,m.^2,-2.*m.*n;
        -m.*n,m.*n,m.^2-n.^2];
      
    % Global transformed reduced stiffness coefficients Q_bar
    Q_bar{i}=inv(T).*Q.*T;
end

Aij=0;
Bij=0;
Dij=0;

for j=1:length(theta)
    A_j=Q_bar{j}*(z(j+1)-z(j));
    Aj{j}=A_j;
    Aij=Aij+Aj{j};
    B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
    Bj{j}=B_j;
    Bij=Bij+Bj{j};
    D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
    Dj{j}=D_j;
    Dij=Dij+Dj{j};
end

A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;

%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];

%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];


%%Shear Extension Coupling
shearf = A\N;
e_xf = shearf(1); e_yf = shearf(2); gamma_xyf = shearf(3);

%%Bending
bendingf = D\M;
k_xf = bendingf(1); k_yf = bendingf(2); k_xyf = bendingf(3);

for k=1:length(theta)
 stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf;%%BOTTOM LAYER
 stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf;%%TOP LAYER
 stress1f(:,k)=stress1; %Bottom
 stress2f(:,k)=stress2; %Top
 b=stress1f(1,:)';
 t=stress2f(1,:)';
 tt=[b t];
 z_1(k,:)=[z(k) z(k+1)];

end
z_1 = z_1';
z_1 = z_1(:);
 figure(1)
 stressx=[stress1f(1,:)' stress2f(1,:)']';
 stressx = stressx(:);
 plot(stressx,z_1)
 xline(0);

 figure(2)
 stressy=[stress1f(2,:)' stress2f(2,:)']';
 stressy = stressy(:);
 plot(stressy,z_1)
 xline(0);

 figure(3)
 stressxy=[stress1f(3,:)' stress2f(3,:)']';
 stressxy = stressxy(:);
 plot(stressxy,z_1)
 xline(0);