i=0;
Pres=[160,200,260,];
x_in_H2=62.5;
x_in_N2=20.8;
x_in_CH4=9.3;
x_in_Ar=3.7;
x_in_NH3=3.7;

N2_in=1;



H2_in=62.5/20.8;
CH4_in=9.3/20.8;
Ar_in=3.7/20.8;
NH3_in=3.7/20.8;
Inertes_in=CH4_in+Ar_in;
Total_in=N2_in+H2_in+NH3_in+Inertes_in;

for cont=1:3
for x=0:0.1:0.9

P=Pres(cont);

N2_out=N2_in-x;
H2_out=H2_in-3*x;
NH3_out=NH3_in+2*x;
Ar_out=Ar_in;
CH4_out=CH4_in;
Inertes_out=Inertes_in;
Total_out=Total_in-2*x;

E1=21500;
E2=47900;
Kp2=sqrt((1/E1)*E2*(1/P^2)*(NH3_out/Total_out)^2/((H2_out/Total_out)^3*(N2_out/Total_out)));
Kp=sqrt((1/P^2)*(NH3_out/Total_out)^2/((H2_out/Total_out)^3*(N2_out/Total_out)));


T=fzero('temp',800,[],[],Kp);
T2=fzero('temp',800,[],[],Kp2);


i=i+1;
Tstore(i)=T-273.14;
Tstore2(i)=T2-273.14;
xstore(i)=x;

end
end

plot(Tstore2(1:10),xstore(1:10),Tstore(1:10),xstore(1:10),Tstore(11:20),xstore(11:20),Tstore(21:30),xstore(21:30))
xlabel('temperatura ºC');
ylabel('conversion %');
legend('P=160,max conv','P=160','P=200','P=260')