I have been looking at this for a couple of days now and have finally convinced myself that this limit doesn't exist and MATLAB is in error. Since you are raising a complex number to a potentially fractional power in this limit expression, and such operations are multi-valued for complex numbers, I would question whether the limit should exist at all simply on that basis. But suppose we restrict ourselves to integer n to avoid that discussion. Start with this expression:

(1 + i/n) ^ (n^2)

Extract the magnitude of the inner part so we are left with the form ( r * u ) ^ (n^2) where u is a unit complex number:

[ sqrt(1 + 1/n^2) * (1 + i/n) / sqrt(1 + 1/n^2) ) ] ^ (n^2)

Separate the magnitude out:

sqrt(1 + 1/n^2) ^ (n^2) * [ (1 + i/n) / sqrt(1 + 1/n^2) ) ] ^ (n^2)

Rewrite the magnitude part using 1/2 exponent instead of sqrt():

[ (1 + 1/n^2) ^ (1/2) ] ^(n^2)

Rearrange:

[ (1 + 1/n^2) ^ (n^2) ] ^(1/2)

The limit of the inside part is clearly e, so we have the magnitude converging to:

e^(1/2)

Now for the unit complex number part, we have a point on the complex unit circle raised to a positive integer power, so the result will clearly also be on the complex unit circle. The question is, does this converge to a particular point or not? To examine that, consider the triangle where it came from and just examine the angle of this triangle:

The angle x associated with our unit complex number is the same as the original triangle it came from, so we have:

tan(x) = (1/n) / 1 = 1/n

So

x = atan(1/n)

The taylor series for atan() is:

taylor(atan(y))

ans =

So for small angles we can just use y, or in our case 1/n.

Since our complex unit vector is raised to the power n^2, we can just multiply that by our small angle to get the expected resulting angle:

(1/n) * n^2 = n

Thus, as n -> infinity, the angle grows without bound.

In conclusion, although the magnitude of the original expression converges to e^(1/2), the angle does not converge, hence the overall limit does not converge. For large n, with a magnitude of e^(1/2) and an angle of n, the original expression (1 + 1i/n)^(n^2) will be pretty close to:

exp(1/2) * (cos(n) + 1i * sin(n)) = exp(1/2 + 1i * n)