diff --git a/src/coordinate_domains.jl b/src/coordinate_domains.jl
index f26e92b..784a886 100644
--- a/src/coordinate_domains.jl
+++ b/src/coordinate_domains.jl
@@ -1,50 +1,45 @@
 """
-    rectangulardomain(Zb, Ze, Δr)
+    rectangulardomain(zbl, ztr, Δr)
 
-Generate initial mesh node coordinates for a rectangular domain ∈ {`Zb`, `Ze`} with
-initial mesh step `Δr`.
-
-`Ze` should be greater than `Zb`.
+Generate initial mesh node coordinates for a rectangular domain extending from the complex
+bottom left corner `zbl` to the top right corner `ztr` with initial mesh step `Δr`.
 """
-function rectangulardomain(Zb::Complex, Ze::Complex, Δr)
-    rZb, iZb = reim(Zb)
-    rZe, iZe = reim(Ze)
+function rectangulardomain(zbl::Complex, ztr::Complex, Δr)
+    rzbl, izbl = reim(zbl)
+    rztr, iztr = reim(ztr)
 
-    rZb > rZe && throw(ArgumentError("real part of Zb greater than Ze"))
-    iZb > iZe && throw(ArgumentError("imag part of Zb greater than Ze"))
+    X = rztr - rzbl
+    Y = iztr - izbl
 
-    X = rZe - rZb
-    Y = iZe - iZb
+    n = ceil(Int, Y/Δr)
+    dy = Y/n
 
-    n = ceil(Int, Y/Δr + 1)
-    dy = Y/(n-1)
-    half_dy = dy/2
+    ## dx = sqrt(Δr² - (dy/2)²), solved for equilateral triangle
+    m = ceil(Int, X/sqrt(Δr^2 - dy^2/4))
+    dx = X/m
+    half_dx = dx/2
 
-    m = ceil(Int, X/sqrt(Δr^2 - dy^2/4) + 1)
-    dx = X/(m-1)
+    T = promote_type(ComplexF64, typeof(zbl), typeof(ztr), typeof(Δr))
+    v = Vector{T}()
+    sizehint!(v, (m+1)*(n+1))
 
-    # precalculate
-    mn = m*n
+    shift = false  # we will displace every other line by dx/2
+    for j = 0:n
+        y = izbl + dy*j
 
-    vlength = mn
-    T = promote_type(typeof(Zb), typeof(Ze), typeof(Δr), Float64)
-    v = Vector{T}()
-    sizehint!(v, vlength)
+        for i = 0:m
+            x = rzbl + dx*i
 
-    on = false
-    @inbounds for j = 0:m-1, i = 0:n-1
-        x = rZb + dx*j
-        y = iZb + dy*i
+            if shift && i == 0
+                continue  # otherwise, we shift out of left bound
+            elseif shift
+                x -= half_dx
+            end
 
-        if (i+1) == n
-            on = !on
-        end
-        if on
-            y += half_dy
-        end
+            if i == m
+                shift = !shift
+            end
 
-        if y < iZe || (abs(y - iZe) < sqrt(eps(real(T))))
-            # Can't just check y <= iZe because of floating point limitations
             push!(v, complex(x, y))
         end
     end
diff --git a/test/DefaultFunction.jl b/test/DefaultFunction.jl
index 56381de..4c03def 100644
--- a/test/DefaultFunction.jl
+++ b/test/DefaultFunction.jl
@@ -15,7 +15,7 @@ xb = -2  # real part begin
 xe = 2  # real part end
 yb = -2  # imag part begin
 ye = 2  # imag part end
-r = 0.2  # initial mesh step
+r = 0.1999  # initial mesh step
 
 origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)
 
@@ -64,11 +64,11 @@ params = GRPFParams(100, maxnodes, 3, maxnodes-1, 1e-9, false)
 @test_logs (:warn,"GRPFParams.maxnodes reached") grpf(defaultfcn, origcoords, PlotData(), params)
 
 # Test with big origcoords
-xb = big"-2"  # real part begin
-xe = big"2"  # real part end
-yb = big"-2"  # imag part begin
-ye = big"2"  # imag part end
-r = big"0.2"  # initial mesh step
+xb = big(xb)  # real part begin
+xe = big(xe)  # real part end
+yb = big(yb)  # imag part begin
+ye = big(ye)  # imag part end
+r = big(r)  # initial mesh step
 
 origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)