I want to evaluate a function defined by the following numerical integral:

```
I1(y_) := NIntegrate((PolyGamma(0, 3*(0.6^2 - x*(1 - x)*y^2)/(4*0.001) + 1)), {x, 0, 1}, PrecisionGoal -> 15, MaxRecursion -> 30)
```

The PolyGamma function is singular if the argument is a negative integer or zero. In this case, singularities only appear after the point $$yapprox1.202$$. You can graphically check it:

```
Manipulate(Plot(PolyGamma(0, 3*(0.6^2 - x*(1 - x)*y^2)/(4*0.001) + 1), {x, 0, 1}), {y, 0, 2})
```

The plot of the integral defined above is:

```
Plot(NIntegrate(Re(PolyGamma(0, 3*(0.6^2 - x*(1 - x)*y^2)/(4*0.001) + 1)), {x, 0, 1}, PrecisionGoal -> 15, MaxRecursion -> 30), {y, 0, 2})
```

If I increase the number of recursions:

```
Plot(NIntegrate(Re(PolyGamma(0, 3*(0.6^2 - x*(1 - x)*y^2)/(4*0.001) + 1)), {x, 0, 1}, PrecisionGoal -> 15, MaxRecursion -> 500), {y, 0, 2})
```

The result is slightly better. The expected plot, roughly speaking, is the “average of these oscilations”(based in physical estimations). I’ve tried to improve the result:

```
Plot(NIntegrate(Re(PolyGamma(0, 3*(0.6^2 - x*(1 - x)*y^2)/(4*0.001) + 1)), {x, 0, 1}, PrecisionGoal -> 15, MaxRecursion -> 500, Method -> "GaussKronrodRule"), {y, 0, 2})
```

But the result is the same. How can I fix the highly oscilatory behaviour? Or at least compute a average of these oscilations.