# sqrt.ode
# newtons method in the complex plane to find square roots of 1
# z^2=1
# x^2-y^2= 1 , 2*x*y=0
# draws  basin boundaries
#
# tolerances
par eps=.001
# discretization of the phase-space
par dx=0,dy=0
# the two roots
par r1=1,s1=0
par r2=-1,s2=0
# the functions
f=x^2-y^2-1
g=2*x*y
# the derivatives of the functions
fx=2*x
fy=-2*y
gx=2*y  
gy=2*x
# the Jacobian -- used in the inverse
det=fx*gy-fy*gx
# a euclidean distance function
dd(x,y)=x*x+y*y
# if close to the root then plot otherwise out of bounds for each root
aux xp[1..2]=if(dd(x-r[j],y-s[j])<eps)then(x0)else(-100)
aux yp[1..2]=if(dd(x-r[j],y-s[j])<eps)then(y0)else(-100)
#
# iteration
x'=x-(f*gy-g*fy)/det
y'=y-(g*fx-f*gx)/det
# initial data
x0'=x0
y0'=y0
# set initial data as parameters
glob 0 t {x0=-1.1+dx*2.2;y0=-1.1+dy*2.2;x=-1.1+dx*2.2;y=-1.1+dy*2.2}
# 
# plot all three sets of points in lousy colors
#
@ meth=discrete,total=20,bound=1000,maxstor=100000,trans=20
@ xp=xp1,yp=yp1,xp2=xp2,yp2=yp2,nplot=2,lt=-1
@ xlo=-1.1,ylo=-1.1,xhi=1.1,yhi=1.1
done

Do a two-parameter range on this ranging over

0 < dx < 1
0 < dy < 1

using the array option