Let us assume that the activation function is a logistic regression denoted as $\sigma()$.
The idea behind cross-entropy (CE) is to optimise the weights $W = [w_1, w_2,...,w_j,...w_k]$ to maximise the log probability - or to minimise the negative log probability.
Here, you are willing to obtain each neuron's derivative of the cost $C^n$ with respect to each of the layers in $W$. Thus, you write $\frac{\partial C}{\partial w_j}$, where $C = [C^1, C^2,...,C^n,...,C^m]$. After some math, which I'll skip here but you can read more about it (in case you're interested here (slide 18 proves useful) and here):
This results in $\frac{1}{n} \sum x_j(\sigma(z)−y)$, where $n$ is the size of your set.
Here, $z=WX+b$, where $X = [x_{11} \ x_{12}...x_{1j}...x_{1k};\quad ....;\quad x_{n1} \ x_{n2}...x_{nj}...x_{nk}]$ ($X$ is an $n$ by $k$ matrix) and $x_{11}..x_{1k}$ are the features you would have per entry, $W$ are the weights as defined above and $b$ is the bias.
In classification, you would like to use this linear dependency of $z$. However, you would want to run it through a non-linear function such as a sigmoid, hereby defined by $\sigma()$ (you can see a proof and read more about it here). $y$ represents the targeted output.
So $w_j$ is the j-th weight of the vector above; $x_j$ is the j-th input vector of an entry, $\sigma(z)$ is the sigmoid applied to the $WX+b$ linear function.
Hope that makes sense.
Note: I am not an expert on backprop, but now having read a bit, I think the following caveat is appropriate. When reading papers or books on neural nets, it is not uncommon for derivatives to be written using a mix of the standard summation/index notation, matrix notation, and multi-index notation (include a hybrid of the last two for tensor-tensor derivatives). Typically the intent is that this should be "understood from context", so you have to be careful!
I noticed a couple of inconsistencies in your derivation. I do not do neural networks really, so the following may be incorrect. However, here is how I would go about the problem.
First, you need to take account of the summation in $E$, and you cannot assume each term only depends on one weight. So taking the gradient of $E$ with respect to component $k$ of $z$, we have
$$E=-\sum_jt_j\log o_j\implies\frac{\partial E}{\partial z_k}=-\sum_jt_j\frac{\partial \log o_j}{\partial z_k}$$
Then, expressing $o_j$ as
$$o_j=\tfrac{1}{\Omega}e^{z_j} \,,\, \Omega=\sum_ie^{z_i} \implies \log o_j=z_j-\log\Omega$$
we have
$$\frac{\partial \log o_j}{\partial z_k}=\delta_{jk}-\frac{1}{\Omega}\frac{\partial\Omega}{\partial z_k}$$
where $\delta_{jk}$ is the Kronecker delta. Then the gradient of the softmax-denominator is
$$\frac{\partial\Omega}{\partial z_k}=\sum_ie^{z_i}\delta_{ik}=e^{z_k}$$
which gives
$$\frac{\partial \log o_j}{\partial z_k}=\delta_{jk}-o_k$$
or, expanding the log
$$\frac{\partial o_j}{\partial z_k}=o_j(\delta_{jk}-o_k)$$
Note that the derivative is with respect to $z_k$, an arbitrary component of $z$, which gives the $\delta_{jk}$ term ($=1$ only when $k=j$).
So the gradient of $E$ with respect to $z$ is then
$$\frac{\partial E}{\partial z_k}=\sum_jt_j(o_k-\delta_{jk})=o_k\left(\sum_jt_j\right)-t_k \implies \frac{\partial E}{\partial z_k}=o_k\tau-t_k$$
where $\tau=\sum_jt_j$ is constant (for a given $t$ vector).
This shows a first difference from your result: the $t_k$ no longer multiplies $o_k$. Note that for the typical case where $t$ is "one-hot" we have $\tau=1$ (as noted in your first link).
A second inconsistency, if I understand correctly, is that the "$o$" that is input to $z$ seems unlikely to be the "$o$" that is output from the softmax. I would think that it makes more sense that this is actually "further back" in network architecture?
Calling this vector $y$, we then have
$$z_k=\sum_iw_{ik}y_i+b_k \implies \frac{\partial z_k}{\partial w_{pq}}=\sum_iy_i\frac{\partial w_{ik}}{\partial w_{pq}}=\sum_iy_i\delta_{ip}\delta_{kq}=\delta_{kq}y_p$$
Finally, to get the gradient of $E$ with respect to the weight-matrix $w$, we use the chain rule
$$\frac{\partial E}{\partial w_{pq}}=\sum_k\frac{\partial E}{\partial z_k}\frac{\partial z_k}{\partial w_{pq}}=\sum_k(o_k\tau-t_k)\delta_{kq}y_p=y_p(o_q\tau-t_q)$$
giving the final expression (assuming a one-hot $t$, i.e. $\tau=1$)
$$\frac{\partial E}{\partial w_{ij}}=y_i(o_j-t_j)$$
where $y$ is the input on the lowest level (of your example).
So this shows a second difference from your result: the "$o_i$" should presumably be from the level below $z$, which I call $y$, rather than the level above $z$ (which is $o$).
Hopefully this helps. Does this result seem more consistent?
Update: In response to a query from the OP in the comments, here is an expansion of the first step.
First, note that the vector chain rule requires summations (see here). Second, to be certain of getting all gradient components, you should always introduce a new subscript letter for the component in the denominator of the partial derivative. So to fully write out the gradient with the full chain rule, we have
$$\frac{\partial E}{\partial w_{pq}}=\sum_i \frac{\partial E}{\partial o_i}\frac{\partial o_i}{\partial w_{pq}}$$
and
$$\frac{\partial o_i}{\partial w_{pq}}=\sum_k \frac{\partial o_i}{\partial z_k}\frac{\partial z_k}{\partial w_{pq}}$$
so
$$\frac{\partial E}{\partial w_{pq}}=\sum_i \left[ \frac{\partial E}{\partial o_i}\left(\sum_k \frac{\partial o_i}{\partial z_k}\frac{\partial z_k}{\partial w_{pq}}\right) \right]$$
In practice the full summations reduce, because you get a lot of $\delta_{ab}$ terms. Although it involves a lot of perhaps "extra" summations and subscripts, using the full chain rule will ensure you always get the correct result.
Best Answer
Firstly, note that $\delta_l$ is nothing else than $\frac{\partial C}{\partial z_l}$, which can be expanded via chain rule to $\frac{\partial C}{\partial a_l}\frac{\partial a_l}{\partial z_l}$.
Moreover, you know that the derivative of each bias can be computed as $\frac{\partial C}{\partial b_l} = \delta_l$ and derivative of each weight is $\frac{\partial C}{\partial w_l} = a_{l-1}\delta_l$. This also holds for the last layer: $\delta_L$ is $\sigma(z)-y$. For derivation of $\delta_L = \frac{\partial C}{\partial z_L} = \sigma(z)-y$, see for example this answer.
You can simply plug it as $\delta_{l+1}$ to compute any preceding $\delta_l$.