Cdd = (7* 927401*10^-29)^2*4*pi*10^-7;
plank = 6626070040 /(10^9*2*pi)* 10^-34;
n = 50000;
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 i=1:4
w = 500*2*pi;
lambda = 0.5*i;
wz=w/(lambda)^(1/3);
wx=wz/lambda;
wy=wx;

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;
end

deltaKONACNO(:,i)=delta;

end

for i=1:4

    plot(N,deltaKONACNO(:,i));
    hold on;

end