For the code shown below there are six values of I for a fixed T. But I want to find out the variation of I (sixth positive value) with T varying from 100 to 500.

gn1=10^(-8);
gp2=10^(-13);
gn2=10^(-7);
gp1=10^(-6);
E1=0.05;
E2=0.1;
N1=10^15;
N2=10^14;
F=2.5*10^19;
Z=10^15;
T=100;
k=8.61*10^(-5);
syms  n p1 p n2 
[n,p1,p,n2]= solve(gn1.*n.*p1+ gp1.*p1*F.*exp(-E1./(k.*T))== gp1.*p.*( N1-p1), gp2.*p.*n2+ gn2.*n2*F.*exp(-E2./(k.*T))== gn2.*n.*( N2-n2), Z== gn1.*n.*p1+ gp2.*p.*n2,((N2-n2)-n)==(( N1-p1)-p))
I= gp2.*p.*n2;