Sorry for the long post. I have tried to reduce the size as much as I can, without the code loosing meaning and would truly appreciate any help. I am essentially running Monte-Carlo Simulations (I have not added the function that actually runs the Monte-Carlo Simulation for Code 1) on the two codes below . Code 2 is just another version of Code 1, but it uses Vector Multiplication in order to increase the speed of the Monte-Carlo Simulation. I am essentially carrying out random sampling from Normal distribution in the functions fibreanglerng(casechoice), vol_fraction(), shearmodulusstrainenergyrng() and diameterng(diamchoice). With Code 2, the Monte Carlo Simulation can be carried out straight away, as the matrices from the aforementioned functions will be configured depending on the number of iterations, whereas Code 1 would require a separate for-loop. Essentially, I am having a problem obtaining the right eCD1 values (I get values which are negative and greater than 1, both of which should not be happening) from Code 2. I have created Code 2 to be almost exactly the same as Code 1 except for the fact that the various Matrix sizes in the code will need to be configured based on the number of iterations (which is where I think I've gone wrong).

Code 1: % Aim of the code is to find eCD1 (which is the equation at the end of the code)
Code:
function [eCD1] = Nielsenneworiginal(casechoice, diamchoice)
symlayup = fibreanglerng(casechoice); % symlaup assigned with 1 x 16 double
opplayup = fliplr(symlayup);
layup = [symlayup,opplayup];

Ef = 235e9; Em = 3.5e9 ; 
vol_f = vol_fraction(); % vol_f assigned with a 1 x 1 double
E11 = (Ef*vol_f) + (Em*(1-vol_f));
E22 = 1/((vol_f/Ef) + ((1-vol_f)/Em));

% Poisson Ratio
nuf = 0.340; num = 0.33;
v12 = (nuf*vol_f) + num*(1-vol_f);

[G12, GIC] = shearmodulusstrainenergyrng(); % G12, GIC assigned with 1 x 1 doubles

T  = 0.000125; % m, Ply thickness

diam = diameterng(diamchoice); % assigned with a 1 x 8 double

NxM = 100e3; % N/m = 10e-6 kN/mm, 100 kN on 100mm = 1kN/mm = 10e6 N/m.
 % Range of 3 standard deviations either side

NyM = 0e3;

NxyM = 0e3;

calcA = 1;% 1 if matrix is calculated
calcB = 0;
calcD = 1;
calcQL = 1;

v21  = (v12*E22)/E11;
Q11  = E11/(1-v12*v21);
Q12  = v21*E11/(1-v12*v21);
Q22  = E22/(1-v12*v21);
Q33  = G12;
Q    = [Q11,Q12,0;Q12,Q22,0;0,0,Q33];

if abs(NxM) < 1e-12
    NxM = 1e-12;
end
if abs(NyM) < 1e-12
    NyM = 1e-12;
end
if abs(NxyM) < 1e-12
    NxyM = 1e-12;
end
W = NxM/NyM;
Y = NxM/NxyM;

E = size(layup,2); % Number of layers
h = E*T; % Thickness of entire laminate
d = 0.25*E; % Number of layers in the above the critical depth
theta=zeros(E);
z=zeros(E+1,1);
x=zeros(d+1,1);

A1=zeros(3,3,E);B1=zeros(3,3,E);D1=zeros(3,3,E);
AS1=zeros(3,3,d);BS1=zeros(3,3,d);DS1=zeros(3,3,d);

%if (z(E+1) ~= (T*E/2)); error('z metric wrong'); end
z(1) = -(T*E)/2;     
for i=1:E
    theta(i)    = (layup(i)*(pi/180));
    z(i+1)      = (z(i)+T);
    m           = (cos(theta(i)));
    n           = (sin(theta(i)));
    Qbar(1,1,i) = (Q(1,1)*(m^4)+2*(Q(1,2)+2*Q(3,3))*(n^2)*(m^2)+Q(2,2)*(n^4));
    Qbar(1,2,i) = ((Q(1,1)+Q(2,2)-4*Q(3,3))*(n^2)*(m^2)+Q(1,2)*((m^4)+(n^4)));
    Qbar(1,3,i) = ((Q(1,1)-Q(1,2)-2*Q(3,3))*(n)*(m^3)+(Q(1,2)-Q(2,2)+2*Q(3,3))*(n^3)*(m));
    Qbar(2,2,i) = (Q(1,1)*(n^4)+2*(Q(1,2)+2*Q(3,3))*(n^2)*(m^2)+Q(2,2)*(m^4));
    Qbar(2,3,i) = ((Q(1,1)-Q(1,2)-2*Q(3,3))*(n^3)*(m)+(Q(1,2)-Q(2,2)+2*Q(3,3))*(n)*(m^3));
    Qbar(3,3,i) = ((Q(1,1)+Q(2,2)-2*Q(1,2)-2*Q(3,3))*(n^2)*(m^2)+Q(3,3)*((m^4)+(n^4)));
    Qbar(2,1,i) = (Qbar(1,2,i));
    Qbar(3,1,i) = (Qbar(1,3,i));
    Qbar(3,2,i) = (Qbar(2,3,i));
    Qbar;
end

A = zeros(3); B = A; D = A; QL = A;
for i= 1:3
    for j= 1:3
        for k= 1:E
            if (calcA); A1(i,j,k) = (Qbar(i,j,k)*(z(k+1)-z(k))); end
            if (calcB); B1(i,j,k) = (1/2)*Qbar(i,j,k)*(z(k+1)^2-z(k)^2); end
            if (calcD); D1(i,j,k) = (1/3)*Qbar(i,j,k)*(z(k+1)^3-z(k)^3); end       
        end
        if (calcA); A(i,j) = sum(A1(i,j,1:E)); end
        if (calcB); B(i,j) = sum(B1(i,j,1:E)); end
        if (calcD); D(i,j) = sum(D1(i,j,1:E)); end
        if (calcQL); QL(i,j) = (sum(Qbar(i,j,1:E)))/E; end
    end
end

A11=A(1,1);A12=A(1,2);A13=A(1,3);A22=A(2,2);A23=A(2,3);A33=A(3,3);A21=A12;...
    A31=A13;A32=A23;
B11=B(1,1);B12=B(1,2);B22=B(2,2);B13=B(1,3);B23=B(2,3);B33=B(3,3);B21=B12;...
    B31=B13;B32=B23;
D11=D(1,1);D12=D(1,2);D13=D(1,3);D22=D(2,2);D23=D(2,3);D33=D(3,3);D21=D12;...
    D31=D13;D32=D23;
QL11=QL(1,1);QL12=QL(1,2);QL13=QL(1,3);QL22=QL(2,2);QL23=QL(2,3);QL33=QL(3,3);...
    QL21=QL12;QL31=QL13;QL32=QL23;

x1(1) = -(T*d)/2;
for i=1:1:0.25*E
    x1(i+1)      = (x1(i)+T);
end

for i= 1:3
    for j= 1:3
        for k= 1:1:d
            if (calcA); AS1(i,j,k) = (Qbar(i,j,k)*(x1(k+1)-x1(k))); end
            if (calcB); BS1(i,j,k) = (1/2)*Qbar(i,j,k)*(x1(k+1)^2-x1(k)^2); end
            if (calcD); DS1(i,j,k) = (1/3)*Qbar(i,j,k)*(x1(k+1)^3-x1(k)^3); end       
        end
        if (calcA); A1S(i,j) = sum(AS1(i,j,1:d)); end
        if (calcB); B1S(i,j) = sum(BS1(i,j,1:d)); end
        if (calcD); D1S(i,j) = sum(DS1(i,j,1:d)); end
    end
end


A1S11=A1S(1,1);A1S12=A1S(1,2);A1S13=A1S(1,3);A1S22=A1S(2,2);A1S23=A1S(2,3);...
    A1S33=A1S(3,3);A1S21=A1S12;A1S31=A1S13;A1S32=A1S23;
B1S11=B1S(1,1);B1S12=B1S(1,2);B1S22=B1S(2,2);B1S13=B1S(1,3);B1S23=B1S(2,3);...
    B1S33=B1S(3,3);B1S21=B1S12;B1S31=B1S13;B1S32=B1S23;
D1S11=D1S(1,1);D1S12=D1S(1,2);D1S13=D1S(1,3);D1S22=D1S(2,2);D1S23=D1S(2,3);...
    D1S33=D1S(3,3);D1S21=D1S12;D1S31=D1S13;D1S32=D1S23;

X1n = (A12^2*A1S13*W + A13^2*A1S12*Y - A12*A13*A1S12*W - A11*A22*A1S13*W +...
    A11*A23*A1S12*W - A12*A23*A1S11*W + A13*A22*A1S11*W - A12*A13*A1S13*Y +...
    A11*A23*A1S13*Y - A13*A23*A1S11*Y - A11*A33*A1S12*Y + A12*A33*A1S11*Y +...
    A23^2*A1S11*W*Y - A12*A23*A1S13*W*Y + A13*A22*A1S13*W*Y - A13*A23*A1S12*W*Y +...
    A12*A33*A1S12*W*Y - A22*A33*A1S11*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...
    A22*A13^2 + A11*A23^2 - A11*A22*A33));
X1d = (A12^2*A1S23*W + A13^2*A1S22*Y - A12*A13*A1S22*W - A12*A23*A1S12*W +...
    A13*A22*A1S12*W - A11*A22*A1S23*W + A11*A23*A1S22*W - A12*A13*A1S23*Y -...
    A13*A23*A1S12*Y + A11*A23*A1S23*Y + A12*A33*A1S12*Y - A11*A33*A1S22*Y +...
    A23^2*A1S12*W*Y - A12*A23*A1S23*W*Y + A13*A22*A1S23*W*Y - A13*A23*A1S22*W*Y +...
    A12*A33*A1S22*W*Y - A22*A33*A1S12*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...
    A22*A13^2 + A11*A23^2 - A11*A22*A33));
Z1n = (A12^2*A1S13*W + A13^2*A1S12*Y - A12*A13*A1S12*W - A11*A22*A1S13*W +...
    A11*A23*A1S12*W - A12*A23*A1S11*W + A13*A22*A1S11*W - A12*A13*A1S13*Y +...
    A11*A23*A1S13*Y - A13*A23*A1S11*Y - A11*A33*A1S12*Y + A12*A33*A1S11*Y +...
    A23^2*A1S11*W*Y - A12*A23*A1S13*W*Y + A13*A22*A1S13*W*Y - A13*A23*A1S12*W*Y +...
    A12*A33*A1S12*W*Y - A22*A33*A1S11*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...
    A22*A13^2 + A11*A23^2 - A11*A22*A33));
Z1d = (A12^2*A1S33*W + A13^2*A1S23*Y - A12*A13*A1S23*W - A12*A23*A1S13*W +...
    A13*A22*A1S13*W + A11*A23*A1S23*W - A11*A22*A1S33*W - A13*A23*A1S13*Y -...
    A12*A13*A1S33*Y + A12*A33*A1S13*Y + A11*A23*A1S33*Y - A11*A33*A1S23*Y +...
    A23^2*A1S13*W*Y - A13*A23*A1S23*W*Y - A12*A23*A1S33*W*Y + A12*A33*A1S23*W*Y +...
    A13*A22*A1S33*W*Y - A22*A33*A1S13*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...
    A22*A13^2 + A11*A23^2 - A11*A22*A33));

if abs(X1n) < 1e-12
    X1n = 1e-12;
end
if abs(X1d) < 1e-12
    X1d = 1e-12;
end
if abs(Z1n) < 1e-12
    Z1n = 1e-12;
end
if abs(Z1d) < 1e-12
    Z1d = 1e-12;
end

X1 = X1n/X1d;

Z1 = Z1n/Z1d;

eCD1 = -(pi^2*(3*D1S11 + 2*D1S12 + 3*D1S22 + 4*D1S33)*(A1S12*A1S23*X1 -...
    A1S13*A1S22*X1 + A1S13*A1S23*Z1 - A1S12*A1S33*Z1 - A1S23^2*X1*Z1 +...
    A1S22*A1S33*X1*Z1))/((0.75*(diam(1))^2)*Z1*(X1 + 1)*(A1S33*A1S12^2 -...
    2*A1S12*A1S13*A1S23 + A1S22*A1S13^2 + A1S11*A1S23^2 - A1S11*A1S22*A1S33));
Code 2:

Code:
function [eCD1] = Nielsennewupdated(casechoice, no_iterations,diamchoice)

[vol_f] = vol_fractionupdated(no_iterations); % vol_f has size no_iterations x 1

symlayup = fibreanglerngupdated(casechoice, no_iterations); 
opplayup = fliplr(symlayup);
layup = [symlayup,opplayup];

Ef = 235e9; Em = 3.5e9 ; %GPa
E11 = (Ef*vol_f) + (Em*(1-vol_f));
E22 = (vol_f/Ef) + ((1-vol_f)/Em);

nuf = 0.340; num = 0.33;
v12 = (nuf*vol_f) + num*(1-vol_f);

% G12, GIC have size no_iterations x 1
[G12, GIC] = shearmodulusstrainenergyrngupdated(no_iterations); 

T  = 0.000125; % m, Ply thickness

% diam has size no_iterations x 8
diam = diameterngupdated(no_iterations,diamchoice);

NxM = 100e3; % N/m = 10e-6 kN/mm, 100 kN on 100mm = 1kN/mm = 10e6 N/m.
 % Range of 3 standard deviations either side

NyM = 0e3;


NxyM = 0e3;


calcA = 1;% 1 if matrix is calculated
calcB = 0;
calcD = 1;
calcQL = 1;

v21  = (v12.*E22)./E11;
Q11  = E11./(1-v12.*v21);
Q12  = v21.*E11./(1-v12.*v21);
Q22  = E22./(1-v12.*v21);
Q33  = G12;
zero_matrix = zeros(no_iterations,1);
Q    = [Q11,Q12,zero_matrix;Q12,Q22,zero_matrix;zero_matrix,zero_matrix,Q33];
Q = Q';


if abs(NxM) < 1e-12
    NxM = 1e-12;
end
if abs(NyM) < 1e-12
    NyM = 1e-12;
end
if abs(NxyM) < 1e-12
    NxyM = 1e-12;
end
W = NxM/NyM;
Y = NxM/NxyM;

E = size(layup,2); % Number of layers
h = E*T; % Thickness of entire laminate
d = 0.25*E; % Number of layers in the above the critical depth
theta=zeros(E);
z=zeros(E+1,1);
x=zeros(d+1,1);

A1=zeros(3,3*no_iterations,E);B1=zeros(3,3*no_iterations,E);... 
    D1=zeros(3,3*no_iterations,E);
A1S=zeros(3,3*no_iterations);B1S=zeros(3,3*no_iterations);...
    D1S=zeros(3,3*no_iterations);      
AS1=zeros(3,3*no_iterations,d);BS1=zeros(3,3*no_iterations,d);...
    DS1=zeros(3,3*no_iterations,d);

%%%%%% Problem could be from here to end of the code %%%%%%%%%%%%%%%%%  


Qbar = zeros(3,3*no_iterations,32);
%if (z(E+1) ~= (T*E/2)); error('z metric wrong'); end
z(1) = -(T*E)/2;     
for i=1:E
    for j=1:3:((3*no_iterations)-2)
    theta(i)    = (layup(i)*(pi/180));
    z(i+1)      = (z(i)+T);
    m           = (cos(theta(i)));
    n           = (sin(theta(i)));
    Qbar(1,j,i) = (Q(1,j)*(m^4)+2*(Q(1,j+1)+2*Q(3,j+2))*(n^2)*(m^2)+Q(2,j+1)*(n^4));
    Qbar(1,j+1,i) = ((Q(1,j)+Q(2,j+1)-4*Q(3,j+2))*(n^2)*(m^2)+Q(1,j+1)*((m^4)+(n^4)));
    Qbar(1,j+2,i) = ((Q(1,j)-Q(1,2)-2*Q(3,3))*(n)*(m^3)+(Q(1,j+1)-Q(2,j+1)+2*Q(3,j+2))*(n^3)*(m));
    Qbar(2,j+1,i) = (Q(1,j)*(n^4)+2*(Q(1,j+1)+2*Q(3,j+2))*(n^2)*(m^2)+Q(2,2)*(m^4));
    Qbar(2,j+2,i) = ((Q(1,j)-Q(1,j+1)-2*Q(3,j+2))*(n^3)*(m)+(Q(1,j+1)-Q(2,j+1)+2*Q(3,j+2))*(n)*(m^3));
    Qbar(3,j+2,i) = ((Q(1,j)+Q(2,j+1)-2*Q(1,j+1)-2*Q(3,j+2))*(n^2)*(m^2)+Q(3,j+2)*((m^4)+(n^4)));
    Qbar(2,j,i) = (Qbar(1,j+1,i));
    Qbar(3,j,i) = (Qbar(1,j+2,i));
    Qbar(3,j+1,i) = (Qbar(2,j+2,i));
    Qbar;
    end
end

A = zeros(3,3*no_iterations); B = A; D = A; QL = A;
for i= 1:3
    for j= 1:3*no_iterations
        for k= 1:E
            if (calcA); A1(i,j,k) = (Qbar(i,j,k)*(z(k+1)-z(k))); 
            end
            if (calcB); B1(i,j,k) = (1/2)*Qbar(i,j,k)*(z(k+1)^2-z(k)^2); end
            if (calcD); D1(i,j,k) = (1/3)*Qbar(i,j,k)*(z(k+1)^3-z(k)^3); end       
        end
        if (calcA); A(i,j) = sum(A1(i,j,1:E)); end
        if (calcB); B(i,j) = sum(B1(i,j,1:E)); end
        if (calcD); D(i,j) = sum(D1(i,j,1:E)); end
        if (calcQL); QL(i,j) = (sum(Qbar(i,j,1:E)))/E; end
    end
end

A11=zeros(1,no_iterations);

i = 0;
for j=1:3:3*no_iterations-2
i = i+1;
    A11(1,i)=A(1,j);A12(1,i)=A(1,j+1);A13(1,i)=A(1,j+2);A22(1,i)=A(2,j+1);...
        A23(1,i)=A(2,j+2);A33(1,i)=A(3,j+2);A21(1,i)=A12(1,i);...
        A31(1,i)=A13(1,i);A32(1,i)=A23(1,i);
    B11(1,i)=B(1,j);B12(1,i)=B(1,j+1);B22(1,i)=B(2,j+1);B13(1,i)=B(1,j+2);...
        B23(1,i)=B(2,j+2);B33(1,i)=B(3,j+2);...
        B21(1,i)=B12(1,i);B31(1,i)=B13(1,i);B32(1,i)=B23(1,i);
    D11(1,i)=D(1,j);D12(1,i)=D(1,j+1);D13(1,i)=D(1,j+2);D22(1,i)=D(2,j+1);...
        D23(1,i)=D(2,j+2);D33(1,i)=D(3,j+2);...
        D21(1,i)=D12(1,i);D31(1,i)=D13(1,i);D32(1,i)=D23(1,i);
    QL11(1,i)=QL(1,j);QL12(1,i)=QL(1,j+1);QL13(1,i)=QL(1,j+2);...
        QL22(1,i)=QL(2,j+1);QL23(1,i)=QL(2,j+2);QL33(1,i)=QL(3,j+2);...
        QL21(1,i)=QL12(1,i);QL31(1,i)=QL13(1,i);QL32(1,i)=QL23(1,i);
end

x1(1) = -(T*d)/2;
for i=1:1:0.25*E
    x1(i+1)      = (x1(i)+T);
end

for i= 1:3
    for j= 1:3*no_iterations
        for k= 1:1:d
            if (calcA); AS1(i,j,k) = (Qbar(i,j,k)*(x1(k+1)-x1(k))); end
            if (calcB); BS1(i,j,k) = (1/2)*Qbar(i,j,k)*(x1(k+1)^2-x1(k)^2); end
            if (calcD); DS1(i,j,k) = (1/3)*Qbar(i,j,k)*(x1(k+1)^3-x1(k)^3); end       
        end
        if (calcA); A1S(i,j) = sum(AS1(i,j,1:d)); end
        if (calcB); B1S(i,j) = sum(BS1(i,j,1:d)); end
        if (calcD); D1S(i,j) = sum(DS1(i,j,1:d)); end
    end
end

i = 0;
for j=1:3:3*no_iterations-2
i = i+1;
    A1S11(1,i)=A1S(1,j);A1S12(1,i)=A1S(1,j+1);A1S13(1,i)=A1S(1,j+2);...
        A1S22(1,i)=A1S(2,j+1);A1S23(1,i)=A1S(2,j+2);...
        A1S33(1,i)=A1S(3,j+2);A1S21(1,i)=A1S12(1,i);A1S31(1,i)=A1S13(1,i);...
        A1S32(1,i)=A1S23(1,i);
    B1S11(1,i)=B1S(1,j);B1S12(1,i)=B1S(1,j+1);B1S22(1,i)=B1S(2,j+1);...
        B1S13(1,i)=B1S(1,j+2);B1S23(1,i)=B1S(2,j+2);...
        B1S33(1,i)=B1S(3,j+2);B1S21(1,i)=B1S12(1,i);B1S31(1,i)=B1S13(1,i);...
        B1S32(1,i)=B1S23(1,i);
    D1S11(1,i)=D1S(1,j);D1S12(1,i)=D1S(1,j+1);D1S13(1,i)=D1S(1,j+2);...
        D1S22(1,i)=D1S(1,j+1);D1S23(1,i)=D1S(2,j+2);...
        D1S33(1,i)=D1S(3,j+2);D1S21(1,i)=D1S12(1,i);D1S31(1,i)=D1S13(1,i);...
        D1S32(1,i)=D1S23(1,i);
end

X1n = (A12(1,:).^2.*A1S13(1,:).*W + A13(1,:).^2.*A1S12(1,:).*Y - ...
    A12(1,:).*A13(1,:).*A1S12(1,:).*W - A11(1,:).*A22(1,:).*A1S13(1,:).*W +...
    A11(1,:).*A23(1,:).*A1S12(1,:).*W - A12(1,:).*A23(1,:).*A1S11(1,:).*W +...
    A13(1,:).*A22(1,:).*A1S11(1,:).*W - A12(1,:).*A13(1,:).*A1S13(1,:).*Y +...
    A11(1,:).*A23(1,:).*A1S13(1,:).*Y - A13(1,:).*A23(1,:).*A1S11(1,:).*Y -...
    A11(1,:).*A33(1,:).*A1S12(1,:).*Y + A12(1,:).*A33(1,:).*A1S11(1,:).*Y +...
    A23(1,:).^2.*A1S11(1,:).*W.*Y - A12(1,:).*A23(1,:).*A1S13(1,:).*W.*Y +...
    A13(1,:).*A22(1,:).*A1S13(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S12(1,:).*W.*Y +...
    A12(1,:).*A33(1,:).*A1S12(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S11(1,:).*W.*Y) ...
    ./(Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) +...
    A22(1,:).*A13(1,:).^2 + A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));
X1d = (A12(1,:).^2.*A1S23(1,:).*W + A13(1,:).^2.*A1S22(1,:).*Y -...
    A12(1,:).*A13(1,:).*A1S22(1,:).*W - A12(1,:).*A23(1,:).*A1S12(1,:).*W +...
    A13(1,:).*A22(1,:).*A1S12(1,:).*W - A11(1,:).*A22(1,:).*A1S23(1,:).*W +...
    A11(1,:).*A23(1,:).*A1S22(1,:).*W - A12(1,:).*A13(1,:).*A1S23(1,:).*Y -...
    A13(1,:).*A23(1,:).*A1S12(1,:).*Y + A11(1,:).*A23(1,:).*A1S23(1,:).*Y +...
    A12(1,:).*A33(1,:).*A1S12(1,:).*Y - A11(1,:).*A33(1,:).*A1S22(1,:).*Y +...
    A23(1,:).^2.*A1S12(1,:).*W.*Y - A12(1,:).*A23(1,:).*A1S23(1,:).*W.*Y +...
    A13(1,:).*A22(1,:).*A1S23(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S22(1,:).*W.*Y +...
    A12(1,:).*A33(1,:).*A1S22(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S12(1,:).*W.*Y)...
    ./(Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) +...
    A22(1,:).*A13(1,:).^2 + A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));
Z1n = (A12(1,:).^2.*A1S13(1,:).*W + A13(1,:).^2.*A1S12(1,:).*Y - A12(1,:).*...
    A13(1,:).*A1S12(1,:).*W - A11(1,:).*A22(1,:).*A1S13(1,:).*W + ...
    A11(1,:).*A23(1,:).*A1S12(1,:).*W - A12(1,:).*A23(1,:).*A1S11(1,:).*W +...
    A13(1,:).*A22(1,:).*A1S11(1,:).*W - A12(1,:).*A13(1,:).*A1S13(1,:).*Y +...
    A11(1,:).*A23(1,:).*A1S13(1,:).*Y - A13(1,:).*A23(1,:).*A1S11(1,:).*Y -...
    A11(1,:).*A33(1,:).*A1S12(1,:).*Y + A12(1,:).*A33(1,:).*A1S11(1,:).*Y +...
    A23(1,:).^2.*A1S11(1,:).*W.*Y - A12(1,:).*A23(1,:).*A1S13(1,:).*W.*Y +...
    A13(1,:).*A22(1,:).*A1S13(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S12(1,:).*W.*Y +...
    A12(1,:).*A33(1,:).*A1S12(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S11(1,:).*W.*Y)./...
    (Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) + A22(1,:).*A13(1,:).^2 +...
    A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));
Z1d = (A12(1,:).^2.*A1S33(1,:).*W + A13(1,:).^2.*A1S23(1,:).*Y -...
    A12(1,:).*A13(1,:).*A1S23(1,:).*W - A12(1,:).*A23(1,:).*A1S13(1,:).*W +...
    A13(1,:).*A22(1,:).*A1S13(1,:).*W + A11(1,:).*A23(1,:).*A1S23(1,:).*W -...
    A11(1,:).*A22(1,:).*A1S33(1,:).*W - A13(1,:).*A23(1,:).*A1S13(1,:).*Y -...
    A12(1,:).*A13(1,:).*A1S33(1,:).*Y + A12(1,:).*A33(1,:).*A1S13(1,:).*Y +...
    A11(1,:).*A23(1,:).*A1S33(1,:).*Y - A11(1,:).*A33(1,:).*A1S23(1,:).*Y +...
    A23(1,:).^2.*A1S13(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S23(1,:).*W.*Y -...
    A12(1,:).*A23(1,:).*A1S33(1,:).*W.*Y + A12(1,:).*A33(1,:).*A1S23(1,:).*W.*Y +...
    A13(1,:).*A22(1,:).*A1S33(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S13(1,:).*W.*Y)./...
    (Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) + A22(1,:).*A13(1,:).^2 +...
    A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));

if abs(X1n) < 1e-12
    X1n = 1e-12;
end
if abs(X1d) < 1e-12
    X1d = 1e-12;
end
if abs(Z1n) < 1e-12
    Z1n = 1e-12;
end
if abs(Z1d) < 1e-12
    Z1d = 1e-12;
end

X1 = X1n/X1d;

Z1 = Z1n/Z1d;

eCD1 = -(pi^2.*(3*D1S11(1,:) + 2.*D1S12(1,:) + 3.*D1S22(1,:) +...
    4.*D1S33(1,:)).*(A1S12(1,:).*A1S23(1,:).*X1 - A1S13(1,:).*A1S22(1,:).*X1...
    +A1S13(1,:).*A1S23(1,:).*Z1 - A1S12(1,:).*A1S33(1,:).*Z1 -...
    A1S23(1,:).^2.*X1.*Z1 + A1S22(1,:).*A1S33(1,:).*X1.*Z1))./...
    ((0.75.*(diam(i,1)).^2).*Z1.*(X1 + 1).*(A1S33(1,:).*A1S12(1,:).^2 -...
    2*A1S12(1,:).*A1S13(1,:).*A1S23(1,:) + A1S22(1,:).*A1S13(1,:).^2 +...
    A1S11(1,:).*A1S23(1,:).^2 - A1S11(1,:).*A1S22(1,:).*A1S33(1,:)))
Thanks a lot