How to compute the cartesian coordinates of each vertices of a Spiral of Theodorus? I'd like to plot integers from 0 to N on a spiral and I have the intuition the easiest way to do this would be to use the Spiral of Theodorus as an approximation.

So, for 0 i would have (0,0), for 1 (1,0), for 2 (1,1) and here begins the difficult part for me.
Can someone help me to figure out a formula or a closure so that I could easily retrieve the coordinates of each positive integer?
 A: Let us start with
$$
(x_k, y_k) = r_k (\cos \phi_k, \sin \phi_k)
$$
you gave 
$$
r_k = \sqrt{k}
$$
so we need $\phi_k$. For the increase we have
$$
\tan \Delta \phi_k = 1 / r_k
$$
and thus
\begin{align}
\phi_1 &= 0 \\
\phi_{k+1} &= \phi_k + \arctan 1/r_k
\end{align}
Here is a construction for your example:

You can fiddle with a live version here: link
A: The coordinates of the $n$th point are
$$\alpha=\sum_{i=1}^n \arctan\left(\frac{1}{\sqrt{i}}\right)$$
$$r=\sqrt{n+1}$$
$$(x,y)=(r\cos(\alpha),r\sin(\alpha))$$
For $n=0$ we get $(1,0)$, for $n=1$ we get $(1,1)$, etc.

A: Approximating $\boldsymbol{\theta_n}$ Using The Euler-Maclaurin Sum Formula
I was not able to compute a closed form for $\theta_n$, but using
$$
\arctan\left(\frac1{\sqrt k}\right)=\frac1{\sqrt{k}}-\frac1{3\sqrt{k}^3}+\frac1{5\sqrt{k}^5}-\frac1{7\sqrt{k}^7}+\frac1{9\sqrt{k}^9}+O\!\left(\frac1{\sqrt{k}^{11}}\right)\tag1
$$
the Euler-Maclaurin Sum Formula says
$$
\begin{align}
\theta_n
&=\sum_{k=2}^n\arctan\left(\frac1{\sqrt k}\right)\\
&=2\sqrt{n}+C+\frac7{6\sqrt{n}}-\frac{41}{120\sqrt{n}^3}+\frac{167}{840\sqrt{n}^5}-\frac{1147}{8064\sqrt{n}^7}+O\!\left(\frac1{\sqrt{n}^9}\right)\tag2
\end{align}
$$
where
$$
C=-2.943181160056894531\tag3
$$
was computed by extending the formula for $\theta_n$ out to $O\!\left(\frac1{\sqrt{n}^{23}}\right)$ and evaluating at $n=100$. The next term in the formula is $\frac{1411}{12672\sqrt{n}^9}\approx\frac1{9\sqrt{n}^9}$ which is a good estimate of the error of the approximation given in $(2)$.
For $n=10$, using the approximation in $(2)$ gives an error in $\theta_n$ of approximately $3.5\times10^{-6}$ radians or $2\times10^{-4}$ degrees. For larger $n$, the error gets smaller.
Then as in the previous answers, we have
$$
(x_n,y_n)=\left(\sqrt{n}\cos(\theta_n),\sqrt{n}\sin(\theta_n)\right)\tag4
$$

Approximating The Curve That Passes Through The Points
Using the first $3$ terms of the approximation, we get
$$
\theta=2r+C+\frac7{6r}\tag5
$$
and solving for $r$, we have
$$
\begin{align}
r
&=\frac{\theta-C+\sqrt{(\theta-C)^2-\frac{28}3}}4\\
&=\frac{\theta-C}2-\frac7{6(\theta-C)}+O\!\left(\frac1{(\theta-C)^3}\right)\tag6
\end{align}
$$
This means that each revolution will increase the radius by $\pi$.
