Python, 143143 134 133 bytes
deff=lambda f(a,b,c,d,R=[-1,1j,1]):
for i in" "*99[R:R=[=[(z:=R[2])-(((a*z+b)*z+c)*z+d)/(z-R[0])/(z-R[1])]+R
for i in" return[x"*99]and[x for x in R[:3]if x.imag==0.0]imag==0]
Attempt This Online!Attempt This Online!
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 to find all approximations, then filtering for whether they're complex.
-9 bytes from @Steffan for various small changes.