Matrix Backpropagation with Softmax and Cross Entropy

I'm having trouble deriving the matrix form of backpropagation. As an example, let's suppose we have the following network: There are two nodes in the input layer plus a bias node fixed at 1, three nodes in the hidden layer plus a bias node fixed at 1, and two output nodes. The signal going into the hidden layer is squashed via the sigmoid function and the signal going into the output layer is squashed via the softmax.

The neural net input and weight matrices would be

$\mathbf{X^1} = \begin{bmatrix} 1 & x_{12} & x_{13} \\ 1 & x_{22} & x_{23} \\ ... & ... & ... \\ 1 & x_{N2} & x_{N3} \end{bmatrix}$

$\mathbf{W^1} = \begin{bmatrix} w_{11} & w_{12} & w_{13} \\ w_{21} & w_{22} & w_{23} \\ w_{31} & w_{32} & w_{33} \end{bmatrix}$

$\mathbf{W^2} = \begin{bmatrix} w_{11} & w_{12} \\ w_{21} & w_{22} \\ w_{31} & w_{32} \\ w_{41} & w_{42} \end{bmatrix}$

and then we can implement the feed-forward process as follows:

$\mathbf{Z^1} = \mathbf{X^1}\mathbf{W^1} = \begin{bmatrix} z_{11} & z_{12} & z_{13} \\ z_{21} & z_{22} & z_{23} \\ ... & ... & ... \\ z_{N1} & z_{N2} & z_{N3} \end{bmatrix}$

$\mathbf{A^1} = sigmoid(Z^1) = \begin{bmatrix} \frac{1}{1 + e^{-z_{11}}} & \frac{1}{1 + e^{-z_{12}}} & \frac{1}{1 + e^{-z_{13}}} \\ \frac{1}{1 + e^{-z_{21}}} & \frac{1}{1 + e^{-z_{22}}} & \frac{1}{1 + e^{-z_{23}}} \\ ... & ... & ... \\ \frac{1}{1 + e^{-z_{N1}}} & \frac{1}{1 + e^{-z_{N2}}} & \frac{1}{1 + e^{-z_{N3}}} \end{bmatrix}$

$\mathbf{X^2} = [\mathbf{1}, \mathbf{A^1}] = \begin{bmatrix} 1 & \frac{1}{1 + e^{-z_{11}}} & \frac{1}{1 + e^{-z_{12}}} & \frac{1}{1 + e^{-z_{13}}} \\ 1 & \frac{1}{1 + e^{-z_{21}}} & \frac{1}{1 + e^{-z_{22}}} & \frac{1}{1 + e^{-z_{23}}} \\ ... & ... & ... & ... \\ 1 & \frac{1}{1 + e^{-z_{N1}}} & \frac{1}{1 + e^{-z_{N2}}} & \frac{1}{1 + e^{-z_{N3}}} \end{bmatrix}$

$\mathbf{Z^2} = \mathbf{X^2}\mathbf{W^2} = \begin{bmatrix} z_{11} & z_{12} \\ z_{21} & z_{22} \\ ... & ... \\ z_{N1} & z_{N2} \end{bmatrix}$

$\mathbf{\hat{Y}} = \mathbf{A^2} = softmax(\mathbf{Z^2}) = \begin{bmatrix} \hat{y}_{11} & \hat{y}_{12} \\ \hat{y}_{21} & \hat{y}_{22} \\ ... & ... \\ \hat{y}_{N1} & \hat{y}_{N2} \end{bmatrix}$

Here, $sigmoid$ is the sigmoid function and similarly $softmax$ is the row-wise softmax function. (Also, note that I've used superscripts to denote the layer of the network for each element.)

My goal is to use gradient descent to find the weight matrices $\mathbf{W^1}$ and $\mathbf{W^2}$ that minimize categorical cross entropy. Can anyone walk me through this process? I have seen it done here, but I can't understand how the solution generalizes to matrix equations.