The problem in your code is
var nextVertex = intersectionPt
which should really be
var nextVertex = [0, 0]
(or at least a constant for all the tiles you generate).
The formula $\sum_{j=0}^4 n_j \mathbf e_j$ for the corner locations produces coordinates in a fixed coordinate system -- it's not relative to the intersection you start out with!
So the gaps you're seeing correspond to the distances between intersection points in the pentagrid, which you're erroneously adding to the tile coordinates.
A different problem that might bother you after you fix that is that the tiles you generate are too large compared to your pentagrid. You seem to be making the spacing between neighboring pentagrind lines equal to side length of the rhombs -- but the average distance between the "ribbons" corresponding to pentagrid lines is clearly longer than that. The tile side length is the width of a ribbon, and there's plenty of space left over between the ribbons!
(Image from https://www.ams.org/publicoutreach/feature-column/fcarc-ribbons)
This is not a problem for generating an infinite tiling, but it will make it difficult to predict which part of the pentagrid to look for intersections in, if you want to find all the tiles in a particular area of the coordinate system.
Curiously, none of your references seem to mention this difference in scale explicitly -- those that give actual formulas seem to assume that the lines on the pentagrid are unit spaced and that the rhomb sides are a unit apart. However, those references also don't actually show the complete pentagrid in the same plane as the finished Penrose tiling. Apparently in their view the pentagrid is only there to define which coordinate sets in $\mathbb Z^5$ correspond to corner points in the tiling, but the actual tiles live in a completely different coordinate system.
If you want each tile corner to appear roughly in the vicinity of the region of the pentagrid that corresponds to it, you need to scale one of the coordinate systems to match. Your sources don't give you the scaling factor to use, so let's derive it. We'll keep the tile side lengths as $1$, and figure out how far apart the pentagrid lines must be.
Suppose we're in the region of the pentagrid in your question with 5-coordinate tuple $(0,0,0,0,0)$, and we now move some large distance $D$ due right. What is the five-tuple corresponding to the region we land in, if the grid spacing is $G$?
Well, we move past $\frac DG$ red lines, crossing them at right angles. The orange lines make an angle of $72^\circ$ with the red ones, so the number of them we cross is an integer close to $\frac DG \cos 72^\circ$. The rounding becomes negligible when we let $D$ got to infinity, so just pretend it is $\frac DG\cos 72^\circ$ exactly. Similarly, the green lines make an angle of $144^\circ$ to the red ones, so the number of them we cross is $\frac DG \cos 144^\circ$. (The factor $\cos 144^\circ$ is negative, corresponding to the fact that we cross these lines in the direction of decreasing region coordinates). And so forth with the cyan and blue, with angles of $216^\circ$ and $288^\circ$.
All in all, we end up with the 5-coordinates, up to rounding
$$ (n_0,n_1,n_2,n_3,n_4) \approx \bigl(\tfrac DG,\tfrac DG\cos 72^\circ,\tfrac DG\cos 144^\circ,\tfrac DG\cos 216^\circ,\tfrac DG\cos 288^\circ\bigr). $$
The coordinates of the tile corner corresponding to that region is
$$ \sum_{j=0}^4 n_j \mathbf e_j = \sum_{j=0}^4 n_j \bigl(\cos(j\cdot 72^\circ), \sin(j\cdot 72^\circ)\bigr). $$
We're just interested in the $x$-coordinate, so plug in the $n_j$s above to get
$$ x \approx \frac DG \sum_{j=0}^4 \cos^2(j\cdot 72^\circ)
= \frac DG \sum_{j=0}^4 \frac{1 - \cos(j\cdot 144^\circ)}2
= \frac DG \cdot \frac{5-\sum_{j=0}^4 \cos(j\cdot 144^\circ)}2
= \frac DG \cdot \frac52 .
$$
(The last equals sign is because the sum of cosines is the change in $x$-coordinate from going all the way around a pentagram with unit sides).
If this $x$ is to grow at the same average rate as $D$, we must set
$$ G = \frac52. $$
Beware, however, that even if you scale the grid by this factor, you still won't get closer than the tiles being "roughly in the vicinity" of the grid intersections that represent them. The nice images in the presentations you link to make it look like each rhombus will contain the grid intersection it comes from -- but it's not possible to align the grid such that this is always the case, nor such that every tile corner is inside the grid region that you compute its location from.
The pentagrid contains places where five grid lines almost but not quite meet in a point, like in this region of your pentagrid:
In this particular almost-intersection the fit is loose enough that we can see it's not exact -- but if your continue the pentagrid far enough, you will eventually find places where almost-intersections like this take place within arbitrarily small areas. The configuration has 6 internal regions and 10 actual intersections, and each of those intersections generates a full-sized tile in the final tiling. The 10 tiles come together in this pattern:
whose diameter is somewhat larger than the pentagrid spacing! Clearly not all of the tiles can lie directly above their defining intersection.
(We can see a differently oriented instance of this pattern in your buggy picture. Those tiles have particularly small gaps between them because the intersections that represent them are close together).
Best Answer
There's a newer proof given in "Two-by-Two Substitution Systems and the Undecidability of the Domino Problem", by Nicolas Ollinger. It seems to be available online at: http://hal.inria.fr/docs/00/26/01/12/PDF/sutica.pdf. Hopefully, it is easier to understand than Robinson's proof.