Actually this has nothing to do with division by zero. When you read the NaN help it lists several ways that NaNs can occur:
"These operations produce NaN:"
- "Any arithmetic operation on a NaN, such as sqrt(NaN)"
- "Addition or subtraction, such as magnitude subtraction of infinities as (+Inf)+(-Inf)"
- "Multiplication, such as 0*Inf"
- "Division, such as 0/0 and Inf/Inf"
- "Remainder, such as rem(x,y) where y is zero or x is infinity"
In your case your code has increasingly large values in both u and v until on the eighth iteration they include several Inf elements (which critically have the same locations in both arrays):
>> any(any(isinf(u)==isinf(v))) % colocated Inf values.
ans = 1
>> any(isnan(u(:)))
ans = 0
>> any(isnan(v(:)))
ans = 0
Thus on the ninth iteration you calculate this: v - ((Delta * u')'), which is thus a subtraction of infinities and (exactly as the documentation states) produces NaNs:
>> tmp = v - (Delta * u')';
>> any(isnan(tmp(:)))
ans = 1