Cdd = (7* 927401*10^-29)^2*4*pi*10^-7;
plank = 6626070040 /(10^9*2*pi)* 10^-34;
n = 50000;
wx = 500* 2 *pi;
wy = 600* 2 *pi;
wz = 700* 2* pi;
M = 166*10^-29 *167;
edd = Cdd/(4*pi)* n^(1/6)*sqrt((M^3*(wx*wy*wz)^(1/3))/plank^5)
eqsystem = @(X)GlobalEquiriblium(X,edd,wx,wy,wz, Cd);

N = linspace(10000,100000,100);
X0 = [1;1;1;90/100;200/100];
for br = 1: 100
    n = N(br);
    edd = Cdd/(4*pi)* n^(1/6)*sqrt((M^3*(wx*wy*wz)^(1/3))/plank^5);
    eqsystem = @(X)GlobalEquiriblium(X,edd,wx,wy,wz,Cd);
    solution = fsolve(eqsystem, X0);
    X0 = solution;
    delta(br) = (solution(5)/solution(4)-1)*100;
    br
end
plot(N,delta);


for i=1:4
w = 500*2*pi;
lambda = 0.5*i;
wz=w/(lambda)^(1/3);
wx=wz/lambda;
wy=wx;
N = linspace(10000,100000,100);
delta = SCN(wx,wy,wz)
plot(N,delta);
hold on
end