# [Python](https://www.python.org), 143 bytes

```python
def f(a,b,c,d,R=[-1,1j,1]):
    for i in" "*99:R=[(z:=R[2])-(((a*z+b)*z+c)*z+d)/(z-R[0])/(z-R[1])]+R
    return[x for x in R[:3]if x.imag==0.0]
````

[Attempt This Online!](https://ato.pxeger.com/run?1=m72soLIkIz9vwYKlpSVpuhY3-1NS0xTSNBJ1knSSdVJ0gmyjdQ11DLN0DGM1rbgUgCAtv0ghUyEzT0lBScvS0gqoQKPKyjYo2ihWU1dDQyNRq0o7SRNIJIOIFE19jSrdoGiDWCgDaEysdhDYoKLUktKivOgKsIkVQBMVgqKtjGMz0xQq9DJzE9NtbQ30DGIhzoK6DuZKAA)

This is the first solution that doesn't use a builtin, and it suffers greatly for it. Instead, I'm using the [Durant-Kerner method](https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method) to find all approximations, then filtering for whether they're complex.