Suppose I have $N$ training examples and there are $K$ classes and the targets have a $1$ of $K$ encoding
Let $t_k^n$ denote the kth component of the nth training target
Let $x^n$ denote the the input to the final hidden layer for the nth training example
Let $a^n = \theta^Tx^n$ denote the activation of the final hidden layer for the nth training example
Let $z^n = \sigma(a^n)$ denote the output of the final hidden layer for the nth training example
Let $y_k^n = \frac{exp(z_k^n)}{\sum_{j=1}^k exp(z_j^n)}$ denote the output of the net for the nth training example
Let $E = -\sum_{n=1}^N \sum_{k=1}^{K} t_k^n \log y_k^n$ denote the loss function
Then: \begin{eqnarray} \frac{\partial}{\partial \theta_{ij}}t_k \log y_k = t_k \frac{d\log y_k}{dy_k} \frac{\partial y_k}{\partial z_j} \frac{dz_j}{da_j} \frac{\partial a_j}{\partial \theta_{ij}} \end{eqnarray} Where:
\begin{eqnarray} \frac{d\log y_k}{dy_k} = \frac{1}{y_k} \end{eqnarray}
\begin{eqnarray} \frac{\partial y_k}{\partial z_j} &=& \frac{\delta_{jk}e^{z_k}(\sum_{l=1}^K e^{z_l}) - e^{z_k}e^{z_j}}{(\sum_{l=1}^K e^{z_l})^2} \\ &=& \frac{\delta_{jk}e^{z_k}}{\sum_{l=1}^K e^{z_l}} - \frac{e^{z_k}}{\sum_{l=1}^K e^{z_l}}\frac{e^{z_j}}{\sum_{l=1}^K e^{z_l}} \\ &=& \delta_{jk}y_k - y_ky_j \\ &=& y_k(\delta_{jk} - y_j) \end{eqnarray}
\begin{eqnarray} \frac{d z_j}{da_j} = \frac{d}{da_j} \sigma(a_j)) = \frac{d}{da_j} \frac{1}{1+e^{-a_j}} = \frac{e^{a_j}}{(1+e^{-a_j})^2} \end{eqnarray}
\begin{eqnarray} \frac{\partial a_j}{\partial \theta_{ij}} = x_i \end{eqnarray}
So:
\begin{eqnarray} \frac{\partial E}{\partial \theta_{ij}} = -\sum_{n=1}^N \sum_{k=1}^{K} t_k^n(\delta_{jk} - y_j^n)\frac{e^{a_j^n}}{(1+e^{-a_j^n})^2}x_i^n \end{eqnarray}
Does this seem right?