Issue with vectorizing code (Matlab)

See more:

MatLab

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)

function [eCD1] = Nielsenneworiginal(casechoice, diamchoice)<br />
symlayup = fibreanglerng(casechoice); % symlaup assigned with 1 x 16 double<br />
opplayup = fliplr(symlayup);<br />
layup = [symlayup,opplayup];<br />
<br />
Ef = 235e9; Em = 3.5e9 ; </pre><br />
vol_f = vol_fraction(); % vol_f assigned with a 1 x 1 double<br />
E11 = (Ef*vol_f) + (Em*(1-vol_f));<br />
E22 = 1/((vol_f/Ef) + ((1-vol_f)/Em));<br />
<br />
% Poisson Ratio<br />
nuf = 0.340; num = 0.33;<br />
v12 = (nuf*vol_f) + num*(1-vol_f);<br />
<br />
[G12, GIC] = shearmodulusstrainenergyrng(); % G12, GIC assigned with 1 x 1 doubles<br />
<br />
T  = 0.000125; % m, Ply thickness<br />
<br />
diam = diameterng(diamchoice); % assigned with a 1 x 8 double<br />
<br />
NxM = 100e3; % N/m = 10e-6 kN/mm, 100 kN on 100mm = 1kN/mm = 10e6 N/m.<br />
 % Range of 3 standard deviations either side<br />
<br />
NyM = 0e3;<br />
<br />
NxyM = 0e3;<br />
<br />
calcA = 1;% 1 if matrix is calculated<br />
calcB = 0;<br />
calcD = 1;<br />
calcQL = 1;<br />
<br />
v21  = (v12*E22)/E11;<br />
Q11  = E11/(1-v12*v21);<br />
Q12  = v21*E11/(1-v12*v21);<br />
Q22  = E22/(1-v12*v21);<br />
Q33  = G12;<br />
Q    = [Q11,Q12,0;Q12,Q22,0;0,0,Q33];<br />
<br />
if abs(NxM) < 1e-12<br />
    NxM = 1e-12;<br />
end<br />
if abs(NyM) < 1e-12<br />
    NyM = 1e-12;<br />
end<br />
if abs(NxyM) < 1e-12<br />
    NxyM = 1e-12;<br />
end<br />
W = NxM/NyM;<br />
Y = NxM/NxyM;<br />
<br />
E = size(layup,2); % Number of layers<br />
h = E*T; % Thickness of entire laminate<br />
d = 0.25*E; % Number of layers in the above the critical depth<br />
theta=zeros(E);<br />
z=zeros(E+1,1);<br />
x=zeros(d+1,1);<br />
<br />
A1=zeros(3,3,E);B1=zeros(3,3,E);D1=zeros(3,3,E);<br />
AS1=zeros(3,3,d);BS1=zeros(3,3,d);DS1=zeros(3,3,d);<br />
<br />
%if (z(E+1) ~= (T*E/2)); error('z metric wrong'); end<br />
z(1) = -(T*E)/2;     <br />
for i=1:E<br />
    theta(i)    = (layup(i)*(pi/180));<br />
    z(i+1)      = (z(i)+T);<br />
    m           = (cos(theta(i)));<br />
    n           = (sin(theta(i)));<br />
    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));<br />
    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)));<br />
    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));<br />
    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));<br />
    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));<br />
    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)));<br />
    Qbar(2,1,i) = (Qbar(1,2,i));<br />
    Qbar(3,1,i) = (Qbar(1,3,i));<br />
    Qbar(3,2,i) = (Qbar(2,3,i));<br />
    Qbar;<br />
end<br />
<br />
A = zeros(3); B = A; D = A; QL = A;<br />
for i= 1:3<br />
    for j= 1:3<br />
        for k= 1:E<br />
            if (calcA); A1(i,j,k) = (Qbar(i,j,k)*(z(k+1)-z(k))); end<br />
            if (calcB); B1(i,j,k) = (1/2)*Qbar(i,j,k)*(z(k+1)^2-z(k)^2); end<br />
            if (calcD); D1(i,j,k) = (1/3)*Qbar(i,j,k)*(z(k+1)^3-z(k)^3); end       <br />
        end<br />
        if (calcA); A(i,j) = sum(A1(i,j,1:E)); end<br />
        if (calcB); B(i,j) = sum(B1(i,j,1:E)); end<br />
        if (calcD); D(i,j) = sum(D1(i,j,1:E)); end<br />
        if (calcQL); QL(i,j) = (sum(Qbar(i,j,1:E)))/E; end<br />
    end<br />
end<br />
<br />
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;...<br />
    A31=A13;A32=A23;<br />
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;...<br />
    B31=B13;B32=B23;<br />
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;...<br />
    D31=D13;D32=D23;<br />
QL11=QL(1,1);QL12=QL(1,2);QL13=QL(1,3);QL22=QL(2,2);QL23=QL(2,3);QL33=QL(3,3);...<br />
    QL21=QL12;QL31=QL13;QL32=QL23;<br />
<br />
x1(1) = -(T*d)/2;<br />
for i=1:1:0.25*E<br />
    x1(i+1)      = (x1(i)+T);<br />
end<br />
<br />
for i= 1:3<br />
    for j= 1:3<br />
        for k= 1:1:d<br />
            if (calcA); AS1(i,j,k) = (Qbar(i,j,k)*(x1(k+1)-x1(k))); end<br />
            if (calcB); BS1(i,j,k) = (1/2)*Qbar(i,j,k)*(x1(k+1)^2-x1(k)^2); end<br />
            if (calcD); DS1(i,j,k) = (1/3)*Qbar(i,j,k)*(x1(k+1)^3-x1(k)^3); end       <br />
        end<br />
        if (calcA); A1S(i,j) = sum(AS1(i,j,1:d)); end<br />
        if (calcB); B1S(i,j) = sum(BS1(i,j,1:d)); end<br />
        if (calcD); D1S(i,j) = sum(DS1(i,j,1:d)); end<br />
    end<br />
end<br />
<br />
<br />
A1S11=A1S(1,1);A1S12=A1S(1,2);A1S13=A1S(1,3);A1S22=A1S(2,2);A1S23=A1S(2,3);...<br />
    A1S33=A1S(3,3);A1S21=A1S12;A1S31=A1S13;A1S32=A1S23;<br />
B1S11=B1S(1,1);B1S12=B1S(1,2);B1S22=B1S(2,2);B1S13=B1S(1,3);B1S23=B1S(2,3);...<br />
    B1S33=B1S(3,3);B1S21=B1S12;B1S31=B1S13;B1S32=B1S23;<br />
D1S11=D1S(1,1);D1S12=D1S(1,2);D1S13=D1S(1,3);D1S22=D1S(2,2);D1S23=D1S(2,3);...<br />
    D1S33=D1S(3,3);D1S21=D1S12;D1S31=D1S13;D1S32=D1S23;<br />
<br />
X1n = (A12^2*A1S13*W + A13^2*A1S12*Y - A12*A13*A1S12*W - A11*A22*A1S13*W +...<br />
    A11*A23*A1S12*W - A12*A23*A1S11*W + A13*A22*A1S11*W - A12*A13*A1S13*Y +...<br />
    A11*A23*A1S13*Y - A13*A23*A1S11*Y - A11*A33*A1S12*Y + A12*A33*A1S11*Y +...<br />
    A23^2*A1S11*W*Y - A12*A23*A1S13*W*Y + A13*A22*A1S13*W*Y - A13*A23*A1S12*W*Y +...<br />
    A12*A33*A1S12*W*Y - A22*A33*A1S11*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...<br />
    A22*A13^2 + A11*A23^2 - A11*A22*A33));<br />
X1d = (A12^2*A1S23*W + A13^2*A1S22*Y - A12*A13*A1S22*W - A12*A23*A1S12*W +...<br />
    A13*A22*A1S12*W - A11*A22*A1S23*W + A11*A23*A1S22*W - A12*A13*A1S23*Y -...<br />
    A13*A23*A1S12*Y + A11*A23*A1S23*Y + A12*A33*A1S12*Y - A11*A33*A1S22*Y +...<br />
    A23^2*A1S12*W*Y - A12*A23*A1S23*W*Y + A13*A22*A1S23*W*Y - A13*A23*A1S22*W*Y +...<br />
    A12*A33*A1S22*W*Y - A22*A33*A1S12*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...<br />
    A22*A13^2 + A11*A23^2 - A11*A22*A33));<br />
Z1n = (A12^2*A1S13*W + A13^2*A1S12*Y - A12*A13*A1S12*W - A11*A22*A1S13*W +...<br />
    A11*A23*A1S12*W - A12*A23*A1S11*W + A13*A22*A1S11*W - A12*A13*A1S13*Y +...<br />
    A11*A23*A1S13*Y - A13*A23*A1S11*Y - A11*A33*A1S12*Y + A12*A33*A1S11*Y +...<br />
    A23^2*A1S11*W*Y - A12*A23*A1S13*W*Y + A13*A22*A1S13*W*Y - A13*A23*A1S12*W*Y +...<br />
    A12*A33*A1S12*W*Y - A22*A33*A1S11*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...<br />
    A22*A13^2 + A11*A23^2 - A11*A22*A33));<br />
Z1d = (A12^2*A1S33*W + A13^2*A1S23*Y - A12*A13*A1S23*W - A12*A23*A1S13*W +...<br />
    A13*A22*A1S13*W + A11*A23*A1S23*W - A11*A22*A1S33*W - A13*A23*A1S13*Y -...<br />
    A12*A13*A1S33*Y + A12*A33*A1S13*Y + A11*A23*A1S33*Y - A11*A33*A1S23*Y +...<br />
    A23^2*A1S13*W*Y - A13*A23*A1S23*W*Y - A12*A23*A1S33*W*Y + A12*A33*A1S23*W*Y +...<br />
    A13*A22*A1S33*W*Y - A22*A33*A1S13*W*Y)/(Y*(A33*A12^2 - 2*A12*A13*A23 +...<br />
    A22*A13^2 + A11*A23^2 - A11*A22*A33));<br />
<br />
if abs(X1n) < 1e-12<br />
    X1n = 1e-12;<br />
end<br />
if abs(X1d) < 1e-12<br />
    X1d = 1e-12;<br />
end<br />
if abs(Z1n) < 1e-12<br />
    Z1n = 1e-12;<br />
end<br />
if abs(Z1d) < 1e-12<br />
    Z1d = 1e-12;<br />
end<br />
<br />
X1 = X1n/X1d;<br />
<br />
Z1 = Z1n/Z1d;<br />
<br />
eCD1 = -(pi^2*(3*D1S11 + 2*D1S12 + 3*D1S22 + 4*D1S33)*(A1S12*A1S23*X1 -...<br />
    A1S13*A1S22*X1 + A1S13*A1S23*Z1 - A1S12*A1S33*Z1 - A1S23^2*X1*Z1 +...<br />
    A1S22*A1S33*X1*Z1))/((0.75*(diam(1))^2)*Z1*(X1 + 1)*(A1S33*A1S12^2 -...<br />
    2*A1S12*A1S13*A1S23 + A1S22*A1S13^2 + A1S11*A1S23^2 - A1S11*A1S22*A1S33)); 

Code 2:

function [eCD1] = Nielsennewupdated(casechoice, no_iterations,diamchoice)<br />
<br />
[vol_f] = vol_fractionupdated(no_iterations); % vol_f has size no_iterations x 1<br />
<br />
symlayup = fibreanglerngupdated(casechoice, no_iterations);<br />
opplayup = fliplr(symlayup);<br />
layup = [symlayup,opplayup];<br />
<br />
Ef = 235e9; Em = 3.5e9 ; %GPa<br />
E11 = (Ef*vol_f) + (Em*(1-vol_f));<br />
E22 = (vol_f/Ef) + ((1-vol_f)/Em);<br />
<br />
nuf = 0.340; num = 0.33;<br />
v12 = (nuf*vol_f) + num*(1-vol_f);<br />
<br />
% G12, GIC have size no_iterations x 1<br />
[G12, GIC] = shearmodulusstrainenergyrngupdated(no_iterations);<br />
<br />
T  = 0.000125; % m, Ply thickness<br />
<br />
% diam has size no_iterations x 8<br />
diam = diameterngupdated(no_iterations,diamchoice);<br />
<br />
NxM = 100e3; % N/m = 10e-6 kN/mm, 100 kN on 100mm = 1kN/mm = 10e6 N/m.<br />
 % Range of 3 standard deviations either side<br />
<br />
NyM = 0e3;<br />
<br />
<br />
NxyM = 0e3;<br />
<br />
<br />
calcA = 1;% 1 if matrix is calculated<br />
calcB = 0;<br />
calcD = 1;<br />
calcQL = 1;<br />
<br />
v21  = (v12.*E22)./E11;<br />
Q11  = E11./(1-v12.*v21);<br />
Q12  = v21.*E11./(1-v12.*v21);<br />
Q22  = E22./(1-v12.*v21);<br />
Q33  = G12;<br />
zero_matrix = zeros(no_iterations,1);<br />
Q    = [Q11,Q12,zero_matrix;Q12,Q22,zero_matrix;zero_matrix,zero_matrix,Q33];<br />
Q = Q&#39;;<br />
<br />
<br />
if abs(NxM) &lt; 1e-12<br />
    NxM = 1e-12;<br />
end<br />
if abs(NyM) &lt; 1e-12<br />
    NyM = 1e-12;<br />
end<br />
if abs(NxyM) &lt; 1e-12<br />
    NxyM = 1e-12;<br />
end<br />
W = NxM/NyM;<br />
Y = NxM/NxyM;<br />
<br />
E = size(layup,2); % Number of layers<br />
h = E*T; % Thickness of entire laminate<br />
d = 0.25*E; % Number of layers in the above the critical depth<br />
theta=zeros(E);<br />
z=zeros(E+1,1);<br />
x=zeros(d+1,1);<br />
<br />
A1=zeros(3,3*no_iterations,E);B1=zeros(3,3*no_iterations,E);...<br />
    D1=zeros(3,3*no_iterations,E);<br />
A1S=zeros(3,3*no_iterations);B1S=zeros(3,3*no_iterations);...<br />
    D1S=zeros(3,3*no_iterations);<br />
AS1=zeros(3,3*no_iterations,d);BS1=zeros(3,3*no_iterations,d);...<br />
    DS1=zeros(3,3*no_iterations,d);<br />
<br />
%%%%%% Problem could be from here to end of the code %%%%%%%%%%%%%%%%%<br />
<br />
<br />
Qbar = zeros(3,3*no_iterations,32);<br />
%if (z(E+1) ~= (T*E/2)); error(&#39;z metric wrong&#39;); end<br />
z(1) = -(T*E)/2;<br />
for i=1:E<br />
    for j=1:3:((3*no_iterations)-2)<br />
    theta(i)    = (layup(i)*(pi/180));<br />
    z(i+1)      = (z(i)+T);<br />
    m           = (cos(theta(i)));<br />
    n           = (sin(theta(i)));<br />
    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));<br />
    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)));<br />
    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));<br />
    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));<br />
    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));<br />
    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)));<br />
    Qbar(2,j,i) = (Qbar(1,j+1,i));<br />
    Qbar(3,j,i) = (Qbar(1,j+2,i));<br />
    Qbar(3,j+1,i) = (Qbar(2,j+2,i));<br />
    Qbar;<br />
    end<br />
end<br />
<br />
A = zeros(3,3*no_iterations); B = A; D = A; QL = A;<br />
for i= 1:3<br />
    for j= 1:3*no_iterations<br />
        for k= 1:E<br />
            if (calcA); A1(i,j,k) = (Qbar(i,j,k)*(z(k+1)-z(k)));<br />
            end<br />
            if (calcB); B1(i,j,k) = (1/2)*Qbar(i,j,k)*(z(k+1)^2-z(k)^2); end<br />
            if (calcD); D1(i,j,k) = (1/3)*Qbar(i,j,k)*(z(k+1)^3-z(k)^3); end<br />
        end<br />
        if (calcA); A(i,j) = sum(A1(i,j,1:E)); end<br />
        if (calcB); B(i,j) = sum(B1(i,j,1:E)); end<br />
        if (calcD); D(i,j) = sum(D1(i,j,1:E)); end<br />
        if (calcQL); QL(i,j) = (sum(Qbar(i,j,1:E)))/E; end<br />
    end<br />
end<br />
<br />
A11=zeros(1,no_iterations);<br />
<br />
i = 0;<br />
for j=1:3:3*no_iterations-2<br />
i = i+1;<br />
    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);...<br />
        A23(1,i)=A(2,j+2);A33(1,i)=A(3,j+2);A21(1,i)=A12(1,i);...<br />
        A31(1,i)=A13(1,i);A32(1,i)=A23(1,i);<br />
    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);...<br />
        B23(1,i)=B(2,j+2);B33(1,i)=B(3,j+2);...<br />
        B21(1,i)=B12(1,i);B31(1,i)=B13(1,i);B32(1,i)=B23(1,i);<br />
    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);...<br />
        D23(1,i)=D(2,j+2);D33(1,i)=D(3,j+2);...<br />
        D21(1,i)=D12(1,i);D31(1,i)=D13(1,i);D32(1,i)=D23(1,i);<br />
    QL11(1,i)=QL(1,j);QL12(1,i)=QL(1,j+1);QL13(1,i)=QL(1,j+2);...<br />
        QL22(1,i)=QL(2,j+1);QL23(1,i)=QL(2,j+2);QL33(1,i)=QL(3,j+2);...<br />
        QL21(1,i)=QL12(1,i);QL31(1,i)=QL13(1,i);QL32(1,i)=QL23(1,i);<br />
end<br />
<br />
x1(1) = -(T*d)/2;<br />
for i=1:1:0.25*E<br />
    x1(i+1)      = (x1(i)+T);<br />
end<br />
<br />
for i= 1:3<br />
    for j= 1:3*no_iterations<br />
        for k= 1:1:d<br />
            if (calcA); AS1(i,j,k) = (Qbar(i,j,k)*(x1(k+1)-x1(k))); end<br />
            if (calcB); BS1(i,j,k) = (1/2)*Qbar(i,j,k)*(x1(k+1)^2-x1(k)^2); end<br />
            if (calcD); DS1(i,j,k) = (1/3)*Qbar(i,j,k)*(x1(k+1)^3-x1(k)^3); end<br />
        end<br />
        if (calcA); A1S(i,j) = sum(AS1(i,j,1:d)); end<br />
        if (calcB); B1S(i,j) = sum(BS1(i,j,1:d)); end<br />
        if (calcD); D1S(i,j) = sum(DS1(i,j,1:d)); end<br />
    end<br />
end<br />
<br />
i = 0;<br />
for j=1:3:3*no_iterations-2<br />
i = i+1;<br />
    A1S11(1,i)=A1S(1,j);A1S12(1,i)=A1S(1,j+1);A1S13(1,i)=A1S(1,j+2);...<br />
        A1S22(1,i)=A1S(2,j+1);A1S23(1,i)=A1S(2,j+2);...<br />
        A1S33(1,i)=A1S(3,j+2);A1S21(1,i)=A1S12(1,i);A1S31(1,i)=A1S13(1,i);...<br />
        A1S32(1,i)=A1S23(1,i);<br />
    B1S11(1,i)=B1S(1,j);B1S12(1,i)=B1S(1,j+1);B1S22(1,i)=B1S(2,j+1);...<br />
        B1S13(1,i)=B1S(1,j+2);B1S23(1,i)=B1S(2,j+2);...<br />
        B1S33(1,i)=B1S(3,j+2);B1S21(1,i)=B1S12(1,i);B1S31(1,i)=B1S13(1,i);...<br />
        B1S32(1,i)=B1S23(1,i);<br />
    D1S11(1,i)=D1S(1,j);D1S12(1,i)=D1S(1,j+1);D1S13(1,i)=D1S(1,j+2);...<br />
        D1S22(1,i)=D1S(1,j+1);D1S23(1,i)=D1S(2,j+2);...<br />
        D1S33(1,i)=D1S(3,j+2);D1S21(1,i)=D1S12(1,i);D1S31(1,i)=D1S13(1,i);...<br />
        D1S32(1,i)=D1S23(1,i);<br />
end<br />
<br />
X1n = (A12(1,:).^2.*A1S13(1,:).*W + A13(1,:).^2.*A1S12(1,:).*Y - ...<br />
    A12(1,:).*A13(1,:).*A1S12(1,:).*W - A11(1,:).*A22(1,:).*A1S13(1,:).*W +...<br />
    A11(1,:).*A23(1,:).*A1S12(1,:).*W - A12(1,:).*A23(1,:).*A1S11(1,:).*W +...<br />
    A13(1,:).*A22(1,:).*A1S11(1,:).*W - A12(1,:).*A13(1,:).*A1S13(1,:).*Y +...<br />
    A11(1,:).*A23(1,:).*A1S13(1,:).*Y - A13(1,:).*A23(1,:).*A1S11(1,:).*Y -...<br />
    A11(1,:).*A33(1,:).*A1S12(1,:).*Y + A12(1,:).*A33(1,:).*A1S11(1,:).*Y +...<br />
    A23(1,:).^2.*A1S11(1,:).*W.*Y - A12(1,:).*A23(1,:).*A1S13(1,:).*W.*Y +...<br />
    A13(1,:).*A22(1,:).*A1S13(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S12(1,:).*W.*Y +...<br />
    A12(1,:).*A33(1,:).*A1S12(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S11(1,:).*W.*Y) ...<br />
    ./(Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) +...<br />
    A22(1,:).*A13(1,:).^2 + A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));<br />
X1d = (A12(1,:).^2.*A1S23(1,:).*W + A13(1,:).^2.*A1S22(1,:).*Y -...<br />
    A12(1,:).*A13(1,:).*A1S22(1,:).*W - A12(1,:).*A23(1,:).*A1S12(1,:).*W +...<br />
    A13(1,:).*A22(1,:).*A1S12(1,:).*W - A11(1,:).*A22(1,:).*A1S23(1,:).*W +...<br />
    A11(1,:).*A23(1,:).*A1S22(1,:).*W - A12(1,:).*A13(1,:).*A1S23(1,:).*Y -...<br />
    A13(1,:).*A23(1,:).*A1S12(1,:).*Y + A11(1,:).*A23(1,:).*A1S23(1,:).*Y +...<br />
    A12(1,:).*A33(1,:).*A1S12(1,:).*Y - A11(1,:).*A33(1,:).*A1S22(1,:).*Y +...<br />
    A23(1,:).^2.*A1S12(1,:).*W.*Y - A12(1,:).*A23(1,:).*A1S23(1,:).*W.*Y +...<br />
    A13(1,:).*A22(1,:).*A1S23(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S22(1,:).*W.*Y +...<br />
    A12(1,:).*A33(1,:).*A1S22(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S12(1,:).*W.*Y)...<br />
    ./(Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) +...<br />
    A22(1,:).*A13(1,:).^2 + A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));<br />
Z1n = (A12(1,:).^2.*A1S13(1,:).*W + A13(1,:).^2.*A1S12(1,:).*Y - A12(1,:).*...<br />
    A13(1,:).*A1S12(1,:).*W - A11(1,:).*A22(1,:).*A1S13(1,:).*W + ...<br />
    A11(1,:).*A23(1,:).*A1S12(1,:).*W - A12(1,:).*A23(1,:).*A1S11(1,:).*W +...<br />
    A13(1,:).*A22(1,:).*A1S11(1,:).*W - A12(1,:).*A13(1,:).*A1S13(1,:).*Y +...<br />
    A11(1,:).*A23(1,:).*A1S13(1,:).*Y - A13(1,:).*A23(1,:).*A1S11(1,:).*Y -...<br />
    A11(1,:).*A33(1,:).*A1S12(1,:).*Y + A12(1,:).*A33(1,:).*A1S11(1,:).*Y +...<br />
    A23(1,:).^2.*A1S11(1,:).*W.*Y - A12(1,:).*A23(1,:).*A1S13(1,:).*W.*Y +...<br />
    A13(1,:).*A22(1,:).*A1S13(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S12(1,:).*W.*Y +...<br />
    A12(1,:).*A33(1,:).*A1S12(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S11(1,:).*W.*Y)./...<br />
    (Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) + A22(1,:).*A13(1,:).^2 +...<br />
    A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));<br />
Z1d = (A12(1,:).^2.*A1S33(1,:).*W + A13(1,:).^2.*A1S23(1,:).*Y -...<br />
    A12(1,:).*A13(1,:).*A1S23(1,:).*W - A12(1,:).*A23(1,:).*A1S13(1,:).*W +...<br />
    A13(1,:).*A22(1,:).*A1S13(1,:).*W + A11(1,:).*A23(1,:).*A1S23(1,:).*W -...<br />
    A11(1,:).*A22(1,:).*A1S33(1,:).*W - A13(1,:).*A23(1,:).*A1S13(1,:).*Y -...<br />
    A12(1,:).*A13(1,:).*A1S33(1,:).*Y + A12(1,:).*A33(1,:).*A1S13(1,:).*Y +...<br />
    A11(1,:).*A23(1,:).*A1S33(1,:).*Y - A11(1,:).*A33(1,:).*A1S23(1,:).*Y +...<br />
    A23(1,:).^2.*A1S13(1,:).*W.*Y - A13(1,:).*A23(1,:).*A1S23(1,:).*W.*Y -...<br />
    A12(1,:).*A23(1,:).*A1S33(1,:).*W.*Y + A12(1,:).*A33(1,:).*A1S23(1,:).*W.*Y +...<br />
    A13(1,:).*A22(1,:).*A1S33(1,:).*W.*Y - A22(1,:).*A33(1,:).*A1S13(1,:).*W.*Y)./...<br />
    (Y.*(A33(1,:).*A12(1,:).^2 - 2*A12(1,:).*A13(1,:).*A23(1,:) + A22(1,:).*A13(1,:).^2 +...<br />
    A11(1,:).*A23(1,:).^2 - A11(1,:).*A22(1,:).*A33(1,:)));<br />
<br />
if abs(X1n) &lt; 1e-12<br />
    X1n = 1e-12;<br />
end<br />
if abs(X1d) &lt; 1e-12<br />
    X1d = 1e-12;<br />
end<br />
if abs(Z1n) &lt; 1e-12<br />
    Z1n = 1e-12;<br />
end<br />
if abs(Z1d) &lt; 1e-12<br />
    Z1d = 1e-12;<br />
end<br />
<br />
X1 = X1n/X1d;<br />
<br />
Z1 = Z1n/Z1d;<br />
<br />
eCD1 = -(pi^2.*(3*D1S11(1,:) + 2.*D1S12(1,:) + 3.*D1S22(1,:) +...<br />
    4.*D1S33(1,:)).*(A1S12(1,:).*A1S23(1,:).*X1 - A1S13(1,:).*A1S22(1,:).*X1...<br />
    +A1S13(1,:).*A1S23(1,:).*Z1 - A1S12(1,:).*A1S33(1,:).*Z1 -...<br />
    A1S23(1,:).^2.*X1.*Z1 + A1S22(1,:).*A1S33(1,:).*X1.*Z1))./...<br />
    ((0.75.*(diam(i,1)).^2).*Z1.*(X1 + 1).*(A1S33(1,:).*A1S12(1,:).^2 -...<br />
    2*A1S12(1,:).*A1S13(1,:).*A1S23(1,:) + A1S22(1,:).*A1S13(1,:).^2 +...<br />
    A1S11(1,:).*A1S23(1,:).^2 - A1S11(1,:).*A1S22(1,:).*A1S33(1,:))); 

Thanks