$ \frac{p_y}{d} = \frac{v_y}{v_z} \implies p_y = \frac{v_y\cdot d}{v_z} = \frac{v_y}{v_z \cdot \tan(\frac{\alpha}{2})}\\ \frac{p_x}{d} = \frac{v_x}{v_x} \implies p_x = \frac{v_x\cdot d}{v_z} = \frac{v_x}{v_x \cdot \tan(\frac{\alpha}{2})} $

$ \frac{p_y}{d} = \frac{v_y}{v_z} \implies p_y = \frac{v_y\cdot d}{v_z} = \frac{v_y}{v_z \cdot \tan(\frac{\alpha}{2})}\\ \frac{p_x}{d} = \frac{v_x}{v_z} \implies p_x = \frac{v_x\cdot d}{v_z} = \frac{v_x}{v_x \cdot \tan(\frac{\alpha}{2})} $

$^1$Assuming you have Column-Major matrices like in OpenGL

Edit

As to how you actually derive it, I'll largely use the explanation byEtay Meiri. It has some additional information and illustration, so you may want to check it, if something seems unclear.

First, your camera is (as mentioned earlier) positioned in the origin. Now assume, you will only project onto a plane of distance $d$ and the plane you project onto is bounded by the aspect ration $ar$ (that is $ar = \frac{window\text{ }width}{window\text{ }height}$ in positive and negative $x$ direction and by $1$ in positive and negative $y$ direction.

You can now determine $d$ by your vertical field of view $\alpha$ (since looking from the side, your camera, the center at the top of your projection plane and the center of the bottom of your projection plane build a triangle and your vertical field of view is the angle at your camera):

$\frac{1}{d} = \tan(\frac{\alpha}{2})\\ \implies d = \frac{1}{\tan(\frac{\alpha}{2})}$

Now you go about calculating projected points. Assume you have an arbitrary point $v = (v_x, v_y, v_z)$, and you want to calculate the point on your plane $p = (p_x, p_y, d)$. Looking at it from the side again, the triangle with your camera, your projection plane top center point and your projected point is similar to the triangle with your camera, your projection plane top center prolonged to have the same z coordinate, as your point $v$ and the point $v$, i.e. they have the same angles$^2$. Therefore, you can now calculate the projected $p_x$ and $p_y$ coordinate:

$ \frac{p_y}{d} = \frac{v_y}{v_z} \implies p_y = \frac{v_y\cdot d}{v_z} = \frac{v_y}{v_z \cdot \tan(\frac{\alpha}{2})}\\ \frac{p_x}{d} = \frac{v_x}{v_x} \implies p_x = \frac{v_x\cdot d}{v_z} = \frac{v_x}{v_x \cdot \tan(\frac{\alpha}{2})} $

To additionally take into account, that the projection plane ranges from $-ar$ to $ar$ in $x$ direction, you would add $ar$ to the denominator in order to make the projected $x$ range be $\left[-1, 1\right]$:

$ p_x = \frac{v_x}{ar \cdot v_x \cdot \tan(\frac{\alpha}{2})} $

If you think about applying this to a matrix now, you need to take into account the division by $z$, which is different for any point with differing $z$ value. Therefore, the division by z is deferred until after the projection matrix is applied (and is done without any code of you, i.e. for OpenGL you'd specify the gl_Position in Clip Space ($\implies$ after projection) and the division is done for you).

The depth test (Z-Test) still needs the $z$ value though, so it needs to be "safed" from the $z$ divide. Therefore you copy your $z$ value to the $w$ component (which is what you got right in your assumption) and end up with the following matrix:

$\left( \begin{array}{cccc} \frac{1}{ar\cdot\tan(\frac{\alpha}{2})}&0&0\\ 0&\frac{1}{\tan(\frac{\alpha}{2})}&0&0\\ 0&0&0&0\\ 0&0&1&0 \end{array} \right)$

The next step is to not project onto a plane, but into a $z$ range of $\left[-1, 1\right]$ (this range is for correct clipping, as is the same range when handling the $x$ and the $y$ coordinate). To be more specific, you don't want all points to end up in that range, but all points between your near and your far plane, so you map $[n, f]$ to $[-1, 1]$. The resulting range is $2$, so you first scale your near-to-far range to $2$. Then you move it to be $[-1, 1]$:

$ f(z) = A\cdot z + B $

And now taking into account that we want to safe the $z$ value from $z$ divide, we get to

$f(z) = A + \frac{B}{z}$

Mapping this scaling to $[-1, 1]$ is a little bit of leg work: you know that any point with $z = n$ (on your near plane) will be projected to $p_z = -1$ and any point with $z = f$ (on your far plane) will be projected to $p_z = 1$. This leads to the following equation system:

$ A + \frac{B}{n} = -1\\ A + \frac{B}{f} = 1 $

Solving this for $A$ and $B$ leads to$^3$:

$ A = \frac{-n-f}{n-f}\\ B = \frac{2fn}{n-f} $

Your third row of the projection matrix must produced the (undivided) $p_z$ (projected) z value. Therefore, you can now choose the individual elements of said row $(\begin{array}{cccc}a&b&c&d\end{array})$ such that the correct $z$ value is produced when multiplying with your point. The correct $p_z$ value is (as established earlier): $p_z = A\cdot z + B$

Therefore this must be the right hand side of your multiplication. Assume your point to project is $(x, y, z, w)$, then you must achieve the following:

$ a \cdot x + b \cdot y + c \cdot z + d \cdot w = A \cdot z + B $

Obviously, your point's $x$ and $y$ coordinate should not influence the projected $z$ coordiante, so you can set $a$ and $b$ to $0$.

$ c \cdot z + d \cdot w = A \cdot z + B $

Since you know that $w$ for any point is $1$, you're left with

$ c \cdot z + d = A \cdot z + B $

And thus, you have $c = A$ and $d = B$.

Now your matrix is complete for projection

$\left( \begin{array}{cccc} \frac{1}{ar\cdot\tan(\frac{\alpha}{2})}&0&0\\ 0&\frac{1}{\tan(\frac{\alpha}{2})}&0&0\\ 0&0&\frac{-n-f}{n-f}&\frac{2fn}{n-f}\\ 0&0&1&0 \end{array} \right)$

Obviously, the tutorial works a little differently in projecting with the vertical field of view, so we can take a look at how you can additionally get to that (with the help of Eric Lengyel:

Instead of projecting onto a plane at $z = d$, we will project onto the near plane, and thus get for a point $v = (v_x, v_y, v_z)$ the projected point $p = (p_x, p_y, n)$:

$ p_x = \frac{v_x \cdot n}{v_z} $

Additionally, you want to map any point with $l\leq x \leq r$ to $\left[-1, 1\right]$, as before. Thus you get

$f(x) = (x-l) \frac{2}{r-l}-1$

Combine those two and you achieve

$p_x = \frac{2n}{r-l}\left(-\frac{v_x}{v_y}\right)- \frac{r+l}{r-l}$

Now put this (and the term for $y$) into the matrix and you get to the first row, first column being $\frac{2n}{r-l}$, and the first row, third column being $\frac{r+l}{r-l}$. Since you may assume $r$ and $l$ to be different, here you can see why those two matrices differ. $r$ and $l$ are guaranteed to be the same in the tutorial's matrix and therefore you get $r+l = ar-ar = 0$ in the denominator, making the third row first column (and second column accordingly) $0$.

$^1$Assuming you have Column-Major matrices like in OpenGL$

$^2$For the $x$ calculation, you will need to take a different point of your projection plane of course, but the idea is the same

$^3$The differences of the sign come from how you orient your camera: in your assumption it is along the negative $z$ axis, whereas in the tutorial it is along the positive $z$ axis