@Carl Love 

Yes, thats exactly it. You can also find all the formulae in the file attached under the original post.

@Carl Love 

Thanks for your reply. That would definitely be an option but my supervidor wants me to use Ritz Minimum Energy Method. The method seems very straight forward until the point of calculating the coefficients. You can see the main porocedure in the follwowing file...

Plate_Buckling.mw

@tomleslie 

G'Day Tom,

Thank you so much for your offered help.

The complete Matlab code is the following. It works fine in Matlab for me.

I really appreciate your efforts.

Best regards

Sam

clear all;

%Materialdata

EL=49627;

ET=15430;

vLT=0.272;

GLT=4800;

G23=4800;

G13=4800;

a=2000;

b=500;

k=5/6;

M=5;

N=5;

R=[-45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45 0 -45 45 0 45 -45]*pi/180;

h=[-24 -23.857 -23.714 -22 -21.857 -21.714 -20 -19.857 -19.714 -18 -17.857 -17.714 -16 -15.857 -15.714 -14 -13.857 -13.714 -12 -11.857 -11.714 -10 -9.857 -9.714 -8 -7.857 -7.714 -6 -5.857 -5.714 -4 -3.857 -3.714 -2 -1.857 -1.714 0 1.714 1.857 2 3.714 3.857 4 5.714 5.857 6 7.714 7.857 8 9.714 9.857 10 11.714 11.857 12 13.714 13.857 14 15.714 15.857 16 17.714 17.857 18 19.714 19.857 20 21.714 21.857 22 23.714 23.857 24];

A(1:3,1:3)=0;

B(1:3,1:3)=0;

D(1:3,1:3)=0;

Q(1:3,1:3)=0;

S(1:3,1:3)=0;

A1(1:2,1:2)=0;

Q1(1:2,1:2)=0;

%Compliance (S-Matrix)

S(1,1)=1/EL;

S(2,2)=1/ET;

S(1,2)=-vLT/EL;

S(3,3)=1/GLT;

S(2,1)=S(1,2);

%Stiffness (Q-matrix) 

Q=inv(S);

Q1(1,1)=G23;

Q1(2,2)=G13;

Q1(1,2)=0;

Q1(2,1)=Q1(1,2);

Q(3,3)=2*Q(3,3);

for i=1:length(R)

T(1,1)=cos(R(i))^2; 

T(2,2)=cos(R(i))^2; 

T(1,2)=sin(R(i))^2; 

T(2,1)=sin(R(i))^2; 

T(3,1)=-sin(R(i))*cos(R(i)); 

T(3,2)=sin(R(i))*cos(R(i)); 

T(1,3)=2*sin(R(i))*cos(R(i));

T(2,3)=-2*sin(R(i))*cos(R(i)); 

T(3,3)=cos(R(i))^2-sin(R(i))^2;

T1(1,1)=cos(R(i)); 

T1(1,2)=-sin(R(i)); 

T1(2,1)=sin(R(i)); 

T1(2,2)=cos(R(i));

Qk(:,:,i)=inv(T)*Q*T; 

Qk(:,3,i)=Qk(:,3,i)/2; 

Qk1(:,:,i)=inv(T1)*Q1*T1;

A=A+Qk(:,:,i)*(h(i+1)-h(i)); 

B=B+.5*Qk(:,:,i)*(h(i+1)^2-h(i)^2); 

D=D+1/3*Qk(:,:,i)*(h(i+1)^3-h(i)^3); 

A1=A1+Qk1(:,:,i)*(h(i+1)-h(i));

end

A; 

B; 

D; 

A1;

matrix1=zeros(M*N,M*N); 

matrix2=zeros(M*N,M*N); 

matrix3=zeros(M*N,M*N); 

matrix4=zeros(M*N,M*N); 

matrix5=zeros(M*N,M*N); 

matrix6=zeros(M*N,M*N); 

matrix7=zeros(M*N,M*N); 

matrix8=zeros(M*N,M*N); 

matrix9=zeros(M*N,M*N);

matrix10=sym(zeros(M*N,M*N)); 

Nxy=sym('Nxy');

for m=1:M

    for n=1:N 

        matrix1(N*(m-1)+n,N*(m-1)+n)=1/4*pi^2*m^2*b/a*D(1,1)+1/4*pi^2*n^2*a/b*D(3,3)+1/4*k*a*b*A1(2,2);

    end

end

matrix1;

for m=1:M

    for n=1:N 

        matrix2(N*(m-1)+n,N*(m-1)+n)=1/4*pi^2*m*n*D(1,2)+1/4*pi^2*m*n*D(3,3);

    end

end

matrix2;

for m=1:M

    for n=1:N 

        matrix3(N*(m-1)+n,N*(m-1)+n)=1/4*k*A1(2,2)*b*m*pi;

    end

end

matrix3;

for m=1:M

    for n=1:N 

    matrix4(N*(m-1)+n,N*(m-1)+n)=1/4*pi^2*m*n*D(1,2)+1/4*pi^2*m*n*D(3,3);

    end

end

matrix4;

for m=1:M

    for n=1:N 

    matrix5(N*(m-1)+n,N*(m-1)+n)=1/4*pi^2*n^2*a/b*D(2,2)+1/4*pi^2*m^2*b/a*D(3,3)+1/4*k*A1(1,1)*a*b;

    end

end

matrix5;

for m=1:M

    for n=1:N matrix6(N*(m-1)+n,N*(m-1)+n)=1/4*k*A1(1,1)*n*pi*a;

    end

end

matrix6;

for m=1:M

    for n=1:N matrix7(N*(m-1)+n,N*(m-1)+n)=1/4*k*A1(2,2)*m*pi*b;

    end

end

matrix7;

for m=1:M

    for n=1:N 

    matrix8(N*(m-1)+n,N*(m-1)+n)=1/4*k*A1(1,1)*n*pi*a;

    end

end

matrix8;

for m=1:M

    for n=1:N matrix9(N*(m-1)+n,N*(m-1)+n)=1/4*k*A1(1,1)*pi^2*n^2*a/b+1/4*k*A1(2,2)*m^2*pi^2*b/a;

    end

end

matrix9;

for m=1:M

    for n=1:N

        for p=1:M

            for q=1:N

                if m==p

                elseif n==q

                elseif rem(m-p,2)==0 

                elseif rem(n-q,2)==0 

                elseif rem(m+p,2)==0 

                elseif rem(n+q,2)==0 

                else

                    i=(m-1)*N+n;

                    j=(p-1)*N+q; 

                    matrix10(i,j)=-8*Nxy*m*n*p*q/((m^2-p^2)*(n^2-q^2));

                end

            end

        end

    end

end

matrix10;

matrixtot=[matrix1 matrix2 matrix3; matrix4 matrix5 matrix6; matrix7 matrix8 matrix9+matrix10];

Nxy=solve(det(matrixtot));

double(Nxy)  

%Computation of Buckling Coefficient

Pansys=35106; % Buckling Coefficient for ANSYS Input 

P=36140; % Buckling Coefficient for Matlab Input

knekkoeff=P*b^2/(pi^2*(D(1,1)*D(2,2))^(1/2)) % Buckling Coefficient 

knekkoeff2=Pansys*b^2/(pi^2*(D(1,1)*D(2,2))^(1/2)) % Buckling Coefficient for ANSYS

DD=(D(1,2)+2*D(3,3))/((D(1,1)*D(2,2))^(1/2))

data=knekkoeff-2*DD % Modified Buckling Coefficient 

data2=knekkoeff2-2*DD % Modified Buckling Coefficient for ANSYS

@Mac Dude 

Thank you very much for your help MacDude. This should get me started.

Regards

Sam