close;
clear;


D_D=0.28;
D_E=0.6;
s1=20;
s1p=0.028;
s2=0.0063;
s3=0.04;
s4=0.8;
s4p=0.027;

D=116.9;
d=1383;
E=3.35;
e=81.6;

q=0;


for count=1:1000,

q=(count-1)/330;

A(1,1)=-s1/(1+s1p*e)-D_D*q^2;
A(1,2)=s2*e;
A(1,3)=0;
A(1,4)=s1*D*s1p/(1+s1p*e)^2+s2*d;
A(2,1)=s1/(1+s1p*e);
A(2,2)=-A(1,2);
A(2,3)=0;
A(2,4)=-A(1,4);
A(3,1)=-s4*e*s4p/(1+s4p*D)^2-s3*E;
A(3,2)=0;
A(3,3)=-s3*D-D_E*q^2;
A(3,4)=s4/(1+s4p*D);
A(4,1)=-A(3,1);
A(4,2)=0;
A(4,3)=s3*D;
A(4,4)=-A(3,4);

Eigvalues=eig(A);
RealE=real(Eigvalues);
MaxE(count)=max(RealE);
qvec(count)=q;

end;

plot(qvec,MaxE,'.r');