It is fairly elementary that close to $x=0$, and in the limit as $\zeta$ becomes large (perhaps hundreds or thousands), $\exp(\zeta x)$ is a good approximation to $(1+x)^\zeta$. But within what range of $x$ is this so? It is easy to find writings on how it is that the exponential function $\exp z$ is
$$\lim_{n→∞}{(1+\frac{z}{n})^n}$$,
but all the texts I have found that expound this seem loathe to give an estimate of just how good an approximation this is, and how small x must be for the approximation to be so-&-so good.
I made an attempt to estimate the error, by considering the expansions for each, looking at which terms in the expansion give the biggest contribution, & by using the identity
$$-(1-x)\ln(1-x)=x(1-\sum_{k=1}^{\infty}\frac{x^k}{k(k+1)})$$
for $|x|<1$ to get the approximation
$$\frac{p^n}{(p)_n}≈\sqrt{1-\frac{n}{p}}\exp(n\sum_{k=1}^\infty{\frac{(\frac{n}{p})^k}{k(k+1)}})$$,
where $(p)_n$ is the descending pochhammer symbol. I'm fairly sure (depending on whether I overlooked any significant terms in usings Stirling's formula (which is alltoo easy to do)), that this expression is pretty good for all n from 0 through p--1. You have to 'cheat' everso slightly at n=p, though, as if you just substitute it the square root bit becomes 0 instead of 1/√(2πp).
The conclusion I came to is that for small $x$, $\exp(\zeta x)$ exceeds $(1+x)^\zeta$ by a factor of approximately
$$\exp(\frac{\zeta x^2}{2})$$.
So if the proportion tolerance be $\epsilon$, then the approximation would be valid in the range bounded by
$$x=±\sqrt{\frac{2\epsilon}{\zeta}}$$,
which I think seems intuitively quite reasonable.
Update
I have been apprised (see comments) of my having used a utterly barbarous method! Overlooking the simple method was a bit barmy, really; but at least there were spin-offs - most particularly that expression for $(p)_n/p^n$. And it shows how with the exponential function | factorial | Pochhammer & all that the threads can just go allover the place & yet always join-up end-to end eventually.
There is a difference between the expression gotten by the simple correct method & the one gotten by my barbarous method. Whereas in the correct method, the series constituting the argument of the exponential is that for
$$x-\ln(1+x)=\sum_{k=2}^\infty\frac{(-x)^k}{k}$$
in my barbarous solution it's
$$x+(1-x)\ln(1-x)=\sum_{k=2}^\infty\frac{x^k}{k(k-1)}$$
The reason for this is when calculating the factor (<1) by which the exponential is multiplied to get the binomial, I just calculated the factor by which the maximum term of the exponential (the term at k (the index of the terms) = $\zeta x$). This gives only the first term - the quadratic - correctly: the series is already wrong at the cubic term, which has a coefficient of +1/6 instead of -1/3. However, the cubic term can be 'restored' by the simple expedient of taking the mean of an exponential weighted by a Gaussian. As the argument ($\zeta x$ in this case) of an exponential becomes large, the distribution of the size of terms with respect to index (k) tends to Gaussian with mean and variance $\zeta x$. So assuming an exponential with stepth $\ln(1-x)$, as in the region of the maximum term, each term is ~(1-x)× the preceding one, we find a correction (and there is no need to do any integration) by noting that the product of the Gaussian & the exponential is
$$\exp(-\frac{k^2}{2\zeta x}+k\ln(1-x))$$
$$=\exp(-\frac{(k-\zeta x\ln(1-x))^2}{2\zeta x}+\frac{\zeta x\ln(1-x)^2}{2})$$
I'll leave out how because we are integrating from -∞ to +∞ with respect to k, it is immaterial that the Gaussian is shifted along the k axis, and we are just left with the correction factor
$$=\exp(\frac{\zeta x\ln(1-x)^2}{2})$$
≈
$$=\exp(\frac{\zeta x^3}{2})$$.
When we divide my barbarous final expression by this (divide rather than multiply, as we have just here been treating the matter reciprocal-wise), the coefficient of the cubic term becomes 1/6 - 1/2 = -1/3 - the correct value!!
I find it interesting that when you use the very crude approximation of taking the peak term of the expansion of the exponential as the typical term, you only get the result correct as far as the quadratic term in the final expression; but that you can extend correctitude as far as the cubic term by the simple expedient of averaging an exponential weighted by a Gaussian.
Welcome to this forum. First of all I want to command you for posting this question, which to me demonstrates that you have a keen interest and enthusiasm for mathematics and the willingness to search for answers by using your own methods. I think it is excellent that you arrive at the correct result, and then manage to rephrase it all (in a humorous way!), in response to my short comment.
The title of your question indicates that you want to explore the relationship between the function $(1+x)^N$ (with $x$ small and $N$ large) and the exponential function $e^{Nx}$.
In the text of the question you chose a slightly better formulation, namely the relation between $(1 + x/N)^N$ and $e^x$, where $x$ is now any real number and $N$ is large. This version is preferable because if we now take $x = -5$, then by taking $N$ large enough we avoid the problem of a negative number raised to a power.
The comparison between these two function is most easily done by expressing both of them as a power of $e$. The first then becomes $e^{N*log(1+x/N)}$. We now write $y = x/N$ and expand the logarithm in a Taylor series. This yields:
$$log(1+y) = y - {y^2}/2 + {y^3}/3 - {y^4}/4 + {y^5}/5...$$
Taking just the first term (which is equivalent to taking the large N limit) yields $e^x$, as expected. Hence the next terms form the correction that we seek. We can write:
$$log(1 + y) - y = -{y^2}/2 * [1 - (2/3)y + (2/4)y^2 - (2/5)y^3 +…]$$
For a crude estimate of the correction it is usually sufficient to consider the first term, which is $-{y^2}/2$. However the simplicity of the series and the fact that successive terms alter in sign, makes it possible to seek even better estimates. If we consider that the factor (2/4) differs only little from (4/9), it is tempting to write the correction term as:
$$log(1+y) - y = -3y^2/(6+4y)$$
This formula works very well. It turns out that the error in this approximation is a roughly a factor $100$ smaller than in the previous one.