Compute backpropagation
You are missing one little step here. At the output node the backpropagation algorithm starts with the loss function $L$. You want to find out how much the weight $w_1$ impacts our loss function thus
$\frac{\partial L}{\partial w_1}$.
We can break this down with the chain rule as
$\frac{\partial L}{\partial w_1} = \frac{\partial L}{\partial f} \frac{\partial f}{\partial w_1} = \frac{\partial L}{\partial f} \frac{\partial f}{\partial \sigma} \frac{\partial \sigma}{\partial w_1}$.
Why do we do this?
We do this to simplify the math. Because the $L$ is a function of $f$ and $f$ is a function of $\sigma$ and $\sigma$ is a function of the weights. But if we ignore the loss function for now.
How do we calculate $\frac{\partial f}{\partial w^2_{i, j}}$?
$\frac{\partial f}{\partial w^2_{i,j}} = \frac{\partial f}{\partial h^2_i} \frac{\partial h^2_i}{\partial w^2_1}$
$f = h_1^2 w_{1,1}^3 + h_2^2 w_{2,1}^3$
$\frac{\partial f}{\partial h^2_i} = w_{i, 1}^3$
$\frac{\partial h^2_i}{\partial w^2_1} = \sigma(W^2h^1) (1 - \sigma(W^2h^1))$
Backpropagation
Based on this answer we can see that $h_i^2$ is a function of the inputs from the previous layer and the activation at the neuron. Thus, we can substitute this function into our equation and continue taking the derivatives and we can see how the error propagates backwards through the network.
$\frac{\partial f}{\partial w^1_{i,j}} = \frac{\partial f}{\partial h^2_i} \frac{\partial h^2_i}{\partial h^1_i} \frac{\partial h^1_i}{\partial w^1_{i,j}}$

(Update)
The answer
This may not be correct, the work needs to be checked.
$f = \begin{bmatrix} w_{11}^3 &amp; w_{21}^3 \end{bmatrix} \begin{bmatrix} h_{1}^2 \\ h_{2}^2 \end{bmatrix} = \textbf{W}^3\textbf{h}^2$
$\textbf{h}^2 = \begin{bmatrix} h_{1}^2 \\ h_{2}^2 \end{bmatrix} = \sigma(\begin{bmatrix} w_{11}^2 &amp; w_{21}^2 \\ w_{21}^2 &amp; w_{22}^2\end{bmatrix} \begin{bmatrix} h_{1}^1 \\ h_{2}^1 \end{bmatrix}) = \sigma(\textbf{W}^2\textbf{h}^1)$
$\textbf{h}^1 = \begin{bmatrix} h_{1}^1 \\ h_{2}^1 \end{bmatrix} = \sigma(\begin{bmatrix} w_{11}^1 &amp; w_{21}^1 \\ w_{21}^1 &amp; w_{22}^1\end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix}) = \sigma(\textbf{W}^1\textbf{X})$
Then using the chain rule
$\frac{\partial f}{\partial \textbf{W}^1} = \frac{\partial f}{\partial \textbf{h}^2} \frac{\partial \textbf{h}^2}{\partial \textbf{h}^1} \frac{\partial \textbf{h}^1}{\partial \textbf{W}^1} = \textbf{W}^3 * \sigma(\textbf{W}^2\textbf{h}^1)(1-\sigma(\textbf{W}^2\textbf{h}^1)) * \sigma(\textbf{W}^1\textbf{X})(1-\sigma(\textbf{W}^1\textbf{X}))$