225.sagews

w, z, r = var('w z r');
a = 1/4;

EE = 1/4*(4-a^2 )*w + w^3 - 1/2 *(2*z - a) * w^2 + a*r - a/2 + 1/4 * a^2 * z - 1/8 * a^3 
EE

-w^2*(2*z - 1/4)/2 + z/64 + w^3 + 63*w/64 + r/4 - 65/512

S = solve(EE, z);
Z = S[0].rhs();
Z

(512*w^3 + 64*w^2 + 504*w + 128*r - 65)/(512*w^2 - 8)

dZ = simplify( diff(Z, w) )
dZ

(1536*w^2 + 128*w + 504)/(512*w^2 - 8) - 1024*w*(512*w^3 + 64*w^2 + 504*w + 128*r - 65)/(512*w^2 - 8)^2

S = solve(dZ==0, w)

Zbp = [ Z(w= S[0].rhs()), Z(w= S[1].rhs()), Z(w= S[2].rhs()), Z(w= S[3].rhs()) ];

N( Zbp[0](r=1) )

-1.75000000000000 - 5.48826389959755e-17*I

rv = 0.00001;
S = [ [N( real(Zbp[0](rv) ) ), N( imag( Zbp[0](rv)) ) ], 
[ N( real(Zbp[1](rv) )), N(imag(Zbp[1](rv) )) ], 
[ N( real(Zbp[2](rv) )), N(imag(Zbp[2](rv) )) ], 
[ N( real(Zbp[3](rv) )), N(imag(Zbp[3](rv) )) ] ];

list_plot( S )

a = animate( [ 
list_plot( 
[ [N( real(Zbp[0](rv) ) ), N( imag( Zbp[0](rv)) ) ], 
[ N( real(Zbp[1](rv) )), N(imag(Zbp[1](rv) )) ], 
[ N( real(Zbp[2](rv) )), N(imag(Zbp[2](rv) )) ], 
[ N( real(Zbp[3](rv) )), N(imag(Zbp[3](rv) )) ] ] )
for rv in srange(0, 1, 0.1) ], 
xmin=-4.1, ymin=-4.1, xmax=4.1, ymax=4.1 )

a.show()