clear all
clc

j = 1;
for k=15:1:50
 syms x;
f = k-16*(0.0625+0.136*x+0.207*x^2+0.278*x^3+0.35*x^4+0.421*x^5+0.493*x^6+0.564*x^7+ 0.635*x^8+0.707*x^9+0.778*x^10+0.85*x^11+0.921*x^12+0.993*x^13+1.064*x^14+1.135*x^15+1.207*x^16+1.278*x^17+1.35*x^18+1.421*x^19+0.00284*25*(1-1*(1-x)^2-2*x*(1-x)^2-3*x^2*(1-x)^2-4*x^3*(1-x)^2-5*x^4*(1-x)^2-6*x^5*(1-x)^2-7*x^6*(1-x)^2-8*x^7*(1-x)^2-9*x^8*(1-x)^2-10*x^9*(1-x)^2-11*x^10*(1-x)^2-12*x^11*(1-x)^2-13*x^12*(1-x)^2-14*x^13*(1-x)^2-15*x^14*(1-x)^2-16*x^15*(1-x)^2-17*x^16*(1-x)^2-18*x^17*(1-x)^2-19*x^18*(1-x)^2-20*x^19*(1-x)^2)/(1-x)^2);
Roots=solve(f);
% Root1 = solve('15-16*(0.0625+0.136*x+0.207*x^2+0.278*x^3+0.35*x^4+0.421*x^5+0.493*x^6+0.564*x^7+ 0.635*x^8+0.707*x^9+0.778*x^10+0.85*x^11+0.921*x^12+0.993*x^13+1.064*x^14+1.135*x^15+1.207*x^16+1.278*x^17+1.35*x^18+1.421*x^19+0.00284*25*(1-1*(1-x)^2-2*x*(1-x)^2-3*x^2*(1-x)^2-4*x^3*(1-x)^2-5*x^4*(1-x)^2-6*x^5*(1-x)^2-7*x^6*(1-x)^2-8*x^7*(1-x)^2-9*x^8*(1-x)^2-10*x^9*(1-x)^2-11*x^10*(1-x)^2-12*x^11*(1-x)^2-13*x^12*(1-x)^2-14*x^13*(1-x)^2-15*x^14*(1-x)^2-16*x^15*(1-x)^2-17*x^16*(1-x)^2-18*x^17*(1-x)^2-19*x^18*(1-x)^2-20*x^19*(1-x)^2)/(1-x)^2)');
%optimal value of LHS=20.6004765804348
% H2/CO=2.3 => LHS=15
% H2/CO=0.5 => LHS=28
rd = double(Roots);
realrd =  sort(rd(find(imag(rd)==0)));
G(1:size(realrd,2),j)=realrd(1:size(realrd,2))';
rroot(j,:)=[k G(1:size(realrd,2),j)'];
%Malpha(j)=1-sqrt(1/k);
j = j+1;
end

plot(1./rroot(:,1),rroot(:,2:end),'-*');
hold on
plot(1./rroot(:,1),ones(size(rroot,1),1),'r:');
hold on
%plot(1./rroot(:,1),Malpha','b:'); 
xlabel('Left hand side value in 6.12','size',14)
ylabel('Value of real roots', 'size',14)
%gtext({'-- Magne alpha2 calculation';'* alpha2 calulation in the appendix'})