p4 = 2;
p5 = 2;
p8 = 1000;
p9 = 12;

syms x2;
x2 = solve(-2-p4+(x2+p4*(x2))/(p9-x2-p4*(x2))*...
    (-1+(p9-x2-p4*(x2))*p8^2/(p8+x2)^2));
x2 = eval(x2);
x1 = (x2+p4*(x2))/p9;
p2 = ((x2+p4.*(x2))./(p9-x2-p4.*(x2)))./exp(x2./(1+x2./p8));

Xfirst = [x2 x1 p2];

%vectors for plot
p3first = [0.13; 0.23];
p2first1 = [p2(1); p2(1)];
p2first2 = [p2(2); p2(2)];
p2first3 = [p2(3); p2(3)];