Another problem with my incrementing

I am having trouble incrementing my dh values all the way through. It's a lot of code, and I've worked on it all, day, but can't find out what I'm not doing right.

%   Fiber
Vf     = 0.6; 
Ef     = 250;               % GPa
rho_f  = 1.8;               % g/(cm^3)
nu_f   = 0.2;
Gf     = Ef / (2*(1 + nu_f)); % GPa
%   Matrix
Vm     = 0.4;
Em     = 3.5;               % GPa
rho_m  = 1.2;               % g/(cm^3)
nu_m   = 0.35;
Gm     = Em / (2*(1 + nu_m)); % GPa
%   Aluminum
E_Al   = 72;   % GPa
rho_Al = 2.78; % g/(cm^3)
nu_Al  = 0.33;
G_Al   = E_Al / (2*(1 + nu_Al)); % GPa
% 3. Please clearly state your assumptions in the development of your
% assessment. Tabulate your material properties in a consistent set of
% units (SI!!).
%
% Part (a)
% Determine the corresponding A, B and D matrix for baseline Al panel and
% each case of composite panels, in which the ply angle in the S2Ep layer
% is assumed to 0, 30, 45, 60, 90 degrees. Each layer thickness =1mm.
% Givens:
%   S2Ep layer
theta = [0 30 45 60 90]; % deg
% GLARE Layout
%   AL
%   S2Ep theta deg.
%   AL
%   S2Ep theta deg.
%   AL
% Calculating Composite (S2Ep) Properties
E_L   = (Vf*Ef) + (Vm*Em);             % GPa
E_T   = (Ef*Em) / ((Ef*Vm) + (Em*Vf)); % GPa
nu_LT = (Vf*nu_f) + (Vm*nu_m);
nu_TL = (E_T*nu_LT) / E_L;
rho_c = (Vf*rho_f) + (Vm*rho_m);
G_LT  = (Gf*Gm) / ((Gf*Vm) + (Gm*Vf)); % GPa
% Calculating Composite (S2Ep) Q Matrix
q11    = E_L / (1-(nu_LT*nu_TL));       % GPa
q22    = E_T / (1-(nu_LT*nu_TL));       % GPa
q12    = nu_LT*E_T / (1-(nu_LT*nu_TL)); % GPa
q66    = G_LT;                          % GPa
Q_S2Ep = [ q11  q12   0;
           q12  q22   0;                % GPa
            0    0   q66];
% Caculating Aluminum Q Matrix
q11_Al = E_Al / (1-(nu_Al^2));          % GPa
q12_Al = nu_Al*q11_Al;                     % GPa
q66_Al = E_Al / (2*(1+nu_Al));          % GPa
Q_Al   = [ q11_Al  q12_Al   0;
           q12_Al  q11_Al   0;                % GPa

Ortho_theta   = [45 -45 -45 45];
  for i = 1:length(Ortho_theta)
     qbar11(i) = (q11.*cosd(Ortho_theta(i)).^4) +...
     (q22.*sind(Ortho_theta(i)).^4) + ...
     (2.*(q12+(2.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2);
     qbar12(i) = (q11 + q22 - ...
     (4.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2 + ...
     (q12.*(cosd(Ortho_theta(i)).^4 + sind(Ortho_theta(i)).^4));
     qbar13(i) = (q11 - q12 - ...
     (2.*q66)).*sind(Ortho_theta(i)).*cosd(Ortho_theta(i)).^3 - ...
     ((q22 - q12 - ...
     (2.*q66)).*(cosd(Ortho_theta(i)).*sind(Ortho_theta(i)).^3));
     qbar21(i) = qbar12(i);
     qbar22(i) = (q11.*sind(Ortho_theta(i)).^4) + ...
     (q22.*cosd(Ortho_theta(i)).^4) + ...
     (2.*(q12 + (2.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2);
      qbar23(i) = (q11 - q12 - ...
     (2.*q66)).*sind(Ortho_theta(i)).^3.*cosd(Ortho_theta(i)) - ...
     ((q22 - q12 - ...
     (2*q66)).*(cosd(Ortho_theta(i)).^3.*sind(Ortho_theta(i))));
      qbar31(i) = qbar13(i);
      qbar32(i) = qbar23(i);
      qbar33(i) = (q11 + q22 - (2*q12) - ...
          (2.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2 + ...
          (q66.*(cosd(Ortho_theta(i)).^4 + sind(Ortho_theta(i)).^4));
      Qbar(1,1,i) = qbar11(i);
      Qbar(1,2,i) = qbar12(i);
      Qbar(1,3,i) = qbar13(i);
      Qbar(2,1,i) = qbar21(i);
      Qbar(2,2,i) = qbar22(i);
      Qbar(2,3,i) = qbar23(i);
      Qbar(3,1,i) = qbar31(i);
      Qbar(3,2,i) = qbar32(i);
      Qbar(3,3,i) = qbar33(i);
  end
  % Calculating Orthotropic A, B, and D Marices
  j = 1;
  dh = [0.01:0.0001:1];
  for i = 1:length(dh)
      Ortho_heights=[-(.7+dh(i)) -.7 0 .7  .7+dh(i)];
      A(:,:,i) = (Ortho_heights(2) - Ortho_heights(1)).*Qbar(:,:,j) + ...
             (Ortho_heights(3) - Ortho_heights(2)).*Qbar(:,:,j+1) + ...
             (Ortho_heights(4) - Ortho_heights(3)).*Qbar(:,:,j+2) + ...
             (Ortho_heights(5) - Ortho_heights(4)).*Qbar(:,:,j+3);
      B(:,:,i) = (1/2)*...
      ((((Ortho_heights(2)).^2 - (Ortho_heights(1)).^2)).*Qbar(:,:,j) + ...
      (((Ortho_heights(3)).^2 - (Ortho_heights(2)).^2)).*Qbar(:,:,j+1) + ...
      (((Ortho_heights(4)).^2 - (Ortho_heights(3)).^2)).*Qbar(:,:,j+2) + ...
      (((Ortho_heights(5)).^2 - (Ortho_heights(4)).^2)).*Qbar(:,:,j+3));
      if B(:,:,i) < 1*10^-6
          B(:,:,i) = 0;
      end
      D(:,:,i) = (1/3)*...
      ((((Ortho_heights(2)).^3 - (Ortho_heights(1)).^3)).*Qbar(:,:,j) + ...
      (((Ortho_heights(3)).^3 - (Ortho_heights(2)).^3)).*Qbar(:,:,j+1) + ...
      (((Ortho_heights(4)).^3 - (Ortho_heights(3)).^3)).*Qbar(:,:,j+2) + ...
      (((Ortho_heights(5)).^3 - (Ortho_heights(4)).^3)).*Qbar(:,:,j+3));
      if ((-(1*10^-6) < D(1,3,i)) && (D(1,3,i) < 1*10^-6))
          flag = 1;
      else
          flag = 0;
      end
      if flag == 1
          breakindex = i;
          break
      end
  end
  i = dh(breakindex)
  fprintf('\nGLARE, [45/-45/-45/45]\n\n');
  fprintf('\nA Matrix\n');
  display(A(:,:,i))
  fprintf('\nB Matrix\n');
  display(B(:,:,i))
  fprintf('\nD Matrix\n');
  display(D(:,:,i))
  ABD_Ortho(:,:,i) = [A(:,:,i) B(:,:,i);
                      B(:,:,i) D(:,:,i)];
  fprintf('------------------------------------------');

[Image Analyst]

Can you be more specific about what the exact problem is? Also, we can't run your code because we don't know what q11, q66, etc. are. Does it run for you? If so, supply enough of the missing code to let it run for us if you want us to help you.

[Thomas]

I'm sorry, I'll update with the q's

[Image Analyst]

And....what about my first question? Are you going to answer that? Or not? If not, I'm done here.

[Thomas]

Basically the iteration starting with the last "for" loop is trying to find the location where D(1,3) and D(2,3) are both equal to, or very close to, 0.

[Geoff Hayes]

Thomas - how does D(2,3) fit into the code as the condition seems only specific to D(1,3,i)

 if ((-(1*10^-6) < D(1,3,i)) && (D(1,3,i) < 1*10^-6))

[Thomas]

I was just testing it for that case first.

[Geoff Hayes]

But where is that line in your code?