Solução iterativa de ponto fixo de uma equação não linear

Em algumas situações nos deparamos com problemas que envolvem a solução de equações não lineares. Por exemplo, qual o valor de $x$ para que $Q(x)$ definido na relação abaixo seja igual a $10^{-3}$?
$$ Q(x) \cong \frac{1}{x\sqrt{2\pi}} (1 - 0,7/x^2) e^{-x^2/2}$$
Até onde eu vejo, uma solução analítica não é possível. Logo, devemos buscar uma solução numérica. Uma forma de resolver esse problema é usar a técnica do ponto fixo. Para isso precisamos encontrar uma expressão do tipo $x_{k+1} = g(x_k)$, com algum valor para $x_0$ que não seja muito distante da solução desejada. Manipulando a expressão $Q(x) = 10^{-3}$ obtemos:
$$ x_{k+1} = \sqrt{-2\ln \left( \frac{10^{-3}x_k\sqrt{2\pi}}{1 - 0,7/x_k^2} \right)}$$
Começando com $x_0 = 2$, obtemos a seguinte sequência
k     x                  Q(x)
 1.   1.1947763   3.1947763
 2.   0.1105417   3.0842345
 3.   0.0096555   3.09389
 4.   0.0008504   3.0930396
 5.   0.0000749   3.0931144
 6.   0.0000066   3.0931078
 7.   0.0000006   3.0931084
 8.   5.104D-08   3.0931084
 9.   4.493D-09   3.0931084

Logo, para $x = 3.0931084$, o valor de $Q(x)$ é igual a $10^{-3}$. Uma observação importante: nem sempre a convergência ocorre e algumas vezes a convergência pode ser muito lenta. Um código Scilab que resolve esse problema é:

clc; 
x = 2;
e = 1;
emax = 1e-8;
k = 0;
pe = 1e-3; p2r = sqrt(2*%pi); pp = pe*p2r;
while e>emax do
    x1 = sqrt(-2*log((pp*x)/(1-0.7/(x*x))));
    e = abs(x - x1);
    x = x1;
    k = k  + 1;
    disp([k, e, x]);
end
qx = exp(-x*x/2)*(1-0.7/(x*x))/(x*p2r);
disp(qx);