Python, 143 134 133 bytes
f=lambda a,b,c,d,R=[-1,1j,1]:[R:=[(z:=R[2])-(((a*z+b)*z+c)*z+d)/(z-R[0])/(z-R[1])]+R for i in" "*99]and[x for x in R[:3]if x.imag==0]
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.