Bell state preparation

As you wrote at point 3, $CNOT$ is the sum of two tensor products, each one involving two matrices. Consider the first tensor product and apply the two matrices respectively to $|0rangle$ and to $|0rangle$ (the $|00rangle$ pair in your state obtained at 2.) and perform the tensor product, so you get the first quarter of the final expansion. Apply the same two matrices respectively to $|0rangle$ and to $|1rangle$ (the $|01rangle$ pair in your state obtained at 2.) and perform the tensor product, so you get the second quarter of the final expansion. Then repeat using the second tensor product of the CNOT expansion, and you get the third and fourth quarters of the final expansion. Shortly, use the distributive property of multiplication.

I would recommend to use a direct matrix representation.

1. An input state $|00rangle$ can be writen as vector $$x= begin{pmatrix} 1 end{pmatrix} $$

2. First step, i.e. Hadamard gate on first qubit and "nothing" on second qubit is described by operation $$ H otimes I = frac{1}{sqrt{2}} begin{pmatrix} 1 & 1 1 & -1 end{pmatrix} otimes begin{pmatrix} 1 & 0 0 & 1 end{pmatrix} = frac{1}{sqrt{2}} begin{pmatrix} 1 & 0 & 1 & 0 0 & 1 & 0 & 1 1 & 0 & -1 & 0 0 & 1 & 0 & -1 end{pmatrix} $$

3. Multiplying vector $x$ with matrix $H otimes I$ gives vector $$ y= frac{1}{sqrt{2}} begin{pmatrix} 1 0 1 end{pmatrix}, $$ which is a first column of the matrix.

4. Now, you can apply CNOT gate described by matrix $$ text{CNOT}= begin{pmatrix} 1 & 0 & 0 & 0 0 & 1 & 0 & 0 0 & 0 & 0 & 1 0 & 0 & 1 & 0 end{pmatrix}, $$ so $$ z=text{CNOT}, y= frac{1}{sqrt{2}} begin{pmatrix} 1 0 0 1 end{pmatrix} $$

In computational basis, vector $z$ is $$ z = frac{1}{sqrt{2}}(|00rangle+|11rangle), $$ which is a desired Bell state.

I will answer the question in a different way. Let's assume your two qubits are represented as below.

Consider your $|00rangle$, represented by $|q_1q_0rangle$, is the state of above circuit before applying any gates, this is same as your original state $|0rangle otimes |0rangle$.

Now applying hadamard in your step 2 is same as applying hadamard on $q_0$. We know $q_0$ is $|0rangle$. So if apply hadamard gate on $|0rangle$ it will be transformed to $frac{1}{sqrt{2}}(|0rangle + |1rangle)$. (check this by doing simple matrix multiplication if you are in doubt).

Let's see $q_1$ now, no gate is applied to it so it remains as $|0rangle$.

So now, the combined state is $frac{1}{sqrt{2}}(|00rangle + |01rangle)$ after applying the hadamard gate. Now apply the CNOT gate (same as your step 3) on the 2 qubits. Classical CNOT gate simply means if control bit is 1 the target bit is flipped, so same logic can be applied for a CNOT gate applied on the qubit. In our case control qubit is $q_0$ and target qubit is $q_1$. So let's apply CNOT on $|q_1q_0rangle$

$CNOT(frac{1}{sqrt{2}}(|00rangle + |01rangle))$ when you apply from rightmost, $|01rangle$ transforms to $|11rangle$ as control qubit ($q_0$) is contributing $|1rangle$ here so the target qubit $q_1$ which is contributing $|0rangle$ flips to $|1rangle$. Similarly if you move left, $|00rangle$ will remain $|00rangle$ as control qubit $q_0$ is contributing $|0rangle$ to this part, so flipping is not required, hence target qubit($q_1$) remains same($|0rangle$). So putting both together you get following

CNOT($frac{1}{sqrt{2}}(|00rangle + |01rangle)$) = $frac{1}{sqrt{2}}(|00rangle + |11rangle)$

This way of thinking would be helpful when you are dealing with bigger transformations and more qubits. Hope it is helpful.

There is another technique to prove the BELL STATE MEASUREMENT.