Well, it is working, but you're operating near a singularity with particularly problematic symmetry.
If you pick a large initial estimate, you're going to end up on the other side of the singularity, with the next step size being even larger. The solution will diverge away from the singularity location. Try a smaller initial estimate, something in the range of (0 0.8) (non-inclusive).
Alternatively, if you want to keep the solver from wandering into bad places, you can try something like this:
f = @(x) 2.5 - 1/x
df = @(x) x.^-2
error = 100
xi = 5 % normally this would cause the solver to diverge
constraint=[0 10]; % but keep the solution constrained
iter = 1
while error > 0.001
xii = xi - (f(xi)/df(xi));
xii = min(max(xii,constraint(1)+eps),constraint(2)-eps)
error = ((xii-xi)/xii)*100;
error = abs(error);
xi=xii;
iter = iter+1;
end
though bear in mind that a crude constraint like this might cause other problems in certain cases.
