I want to perform the following integration: $ $ f(x)=\int_{0.01}^{x-0.01} \int_{0.1}^{10} \int_{0.1}^{10} \frac{8(x-x_1)^2Ep^2}{2x_1^2 Ep^2 s_2^2+2(x-x_1)^2 Ep^2 s_1s_2-s_1s_2^2-4x_1(x-x_1)Ep^2 \sqrt{s_1}\sqrt{s_2}s_2}ds_2ds_1dx_1,$ $ and then plot $ f(x)$ for $ x \in[0,1]$ . In the above expression, $ Ep$ is given. Such computation cannot be performed analytically, thus one has to take the numerical approach.

I wrote the following lines of codes in Mathematica but apparently do not work.

Ep=10; f[x_, x1_, s1_, s2_] := (8*(x - x1)^2*Ep^2)/(2*(x1^2*Ep^2*s2^2) + 2*((x - x1)^2*Ep^2*s1*s2) - (s1*s2^2) -4*x1*(x - x1)*(Ep^2*Sqrt[s1]*Sqrt[s2]*s2)) Ng[x_?NumericQ] := NIntegrate[f[x, x1, s1, s2], {x1, 0.01, x-0.01}, {s1, 0.1, 10}, {s2, 0.1, 10}] Plot[Ng[x], {x,0.001, 1}] 

Does anyone know a clever way to do this integration?