Ex 2.4

clear
syms rs;
xtp = 0.5 + pi*rs^2*(0.5 - rs)/(1 - pi*rs^2);
xtpprim = diff(xtp,rs);
sol = solve(xtpprim,rs);
nsol = eval(sol);

c = [pi,0,-3,1];
r = roots(c);

% hitta rätt rot jämför resultaten
nsol1 = nsol(nsol>0 & nsol<0.5);
r1 = r(r>0 & r<0.5);
skillnad = abs(nsol1 - r1);
disp(['Skillnaden blev ',num2str(skillnad)]);

% Räkna fram xtp och skriv ut resultaten
xtp1 = subs(xtp,rs,nsol1);
disp(['xtp=',num2str(xtp1)]);
disp(['r=',num2str(r1)]);
Skillnaden blev 1.1102e-16
xtp=0.60107
r=0.40071