Received this error, not sure how to fix

**This code was copied from a textbook so I'm not sure why it does not work.

Index exceeds the number of array elements. Index must not exceed 6.

Error in HW6 (line 62)

wmax(i,j)=wmax(i,j)+(16*PO/(pi^6))*sin(pi*M(mi)*xc/a)*sin(pi*N(ni)*yc/b)/(M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4));

Here is the code:

E1=38.6E9;
E2=8.27E9;
G12=4.14E9;
v12=0.26;
v21=E2*v12/E1;
theta=[0 90 90 0];
h=[-0.0025 -0.00125 0 0.00125 0.0025];
a=0.5;
b=0.25;
M=1:2:11;
N=1:2:11;
PO=10;
xc=a/2;
yc=b/2;
z=h(1);
x=0:0.01:xc;
D=zeros(3,3);
wmax=zeros(6,6);
WO=zeros(1,length(x));
epsilonX=zeros(1,length(x));
epsilonY=zeros(1,length(x));
epsilonXY=zeros(1,length(x));
for i=1:length(theta)
    m=cosd(theta(i));
    n=sind(theta(i));
    Q11=E1/(1-v12*v21);
    Q12=E2*v12/(1-v12*v21);
    Q22=E2/(1-v12*v21);
    Q66=G12;
    T1=[m^2 n^2 2*m*n;n^2 m^2 -2*m*n;-m*n m*n m^2-n^2];
    T2=[m^2 n^2 m*n;n^2 m^2 -m*n;-2*m*n 2*m*n m^2-n^2];
    Q=[Q11 Q12 0;Q12 Q22 0;0 0 Q66];
    QBAR(:,:,i)=inv(T1)*Q*T2;
    Di(:,:,i)=(1/3)*(h(i+1)^3-h(i)^3)*QBAR(:,:,i);
end
for j=1:length(theta)
    D=Di(:,:,j)+D;
end
D11=D(1,1);
D12=D(1,2);
D22=D(2,2);
D66=D(3,3);
for i=1:length(M)
    for j=1:length(N)
        for mi=1:i
            for ni=1+j
                wmax(i,j)=wmax(i,j)+(16*PO/(pi^6))*sin(pi*M(mi)*xc/a)*sin(pi*N(ni)*yc/b)/(M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4));
            end
        end
    end
end

Index exceeds the number of array elements. Index must not exceed 6.

QBARTOP=QBAR(:,:,1);
for i=1:length(x)
    for mi=1:length(M)
        for ni=1:length(N)
            WO(i)=WO(i)+(16*PO/(pi^6))*sin(pi*M(mi)*x(i)/a)*sin(pi*N(ni)*yc/b)/(M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4));
            
            epsilonX(i)=epsilonX(i)-z*(-(16*PO/(pi^6))*(pi*M(mi)/a)^2*sin(pi*M(mi)*x(i)/a)*sin(pi*N(ni)*yc/b)/(M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4)));
            epsilonY(i)=epsilonY(i)-z*(-(16*PO/(pi^6))*(pi*N(ni)/b)^2*sin(pi*M(mi)*x(i)/a)*sin(pi*N(ni)*yc/b)/(M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4)));
            epsilonXY(i)=epsilonXY(i)-2*z*((16*PO/(pi^6))*(pi*N(ni)/b)*(pi*M(mi)/a)*cos(pi*M(mi)*x(i)/a)*cos(pi*N(ni)*yc/b)/M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4));
            deltaX(i)=QBARTOP(1,1)*epsilonX(i)+QBARTOP(1,2)*epsilonY(i)+QBARTOP(1,3)*epsilonXY(i);
            deltaY(i)=QBARTOP(1,2)*epsilonX(i)+QBARTOP(2,2)*epsilonY(i)+QBARTOP(2,3)*epsilonXY(i);
        end
    end
end
subplot(2,1,1);
plot(x,WO, 'r.-'); grid on;
xlabel('Distance from edge x, m');
ylabel('Deflection wo, m');
subplot(2,1,2);
plot(x,abs(deltaX), 'g-*') ; grid on; hold on;
plot(x,abs(deltaY),'k-o');
xlabel('Distance from edge x, m');
ylabel("Bending stress, N/m^2");

Answer 1

[Torsten]

Maybe you mean

for ni=1:j

instead of

for ni=1+j
    

If not, you try to access N(7) which does not exist.

[Anita Burns]

Thank you!! I have looked over this so many times and not caught that. I feel silly.

Answer 2

[Walter Roberson]

M=1:2:11;
N=1:2:11;

those are length 6

for i=1:length(M)
    for j=1:length(N)
        for mi=1:i
            for ni=1+j

when j becomes length(N) then ni becomes 1+length(N)

                wmax(i,j)=wmax(i,j)+(16*PO/(pi^6))*sin(pi*M(mi)*xc/a)*sin(pi*N(ni)*yc/b)/(M(mi)*N(ni)*(D11*(M(mi)/a)^4+2*(D12+2*D66)*(M(mi)/a)^2*(N(ni)/b)^2+D22*(N(ni)/b)^4));

You index N(ni) when ni is potentially 1+length(N) . So you are indexing N(1+length(N)) which is going to be out of range.