How does the slow All-pairs-shortest-paths algorithm work?
Let me tell you what Bellman-Fords is actually doing because the sentence by relaxing all edges once for every vertex in the graph is not very accurate. In Bellman-Ford you care about single-source-shortest-path. Since every path in the graph have at most $n$ nodes it has at most $n-1$ edges. 

Let's say the source node is $s$, then create a table $distance[\cdot]$ such that $distance[s] = 0$ and it is $\infty$ elsewhere. Let's solve the following dynamic programming problem, compute shortest path from $s$ to $u$ (for every $u$) that uses at at least $t$ edges.

The table we already filled is correct for the case of $t=0$. If you have this table already computed for some $t$ after relaxing all edges (in any order) it will be correct for $t+1$, since you extend all shortest path that uses $t$ edges with at least one more edge by relaxing the end of the path. Notice that you don't need to relax every edge $n$ times, but $n - 1$ times since largest path will have at most $n - 1$ edges.

The algorithm you showed is doing that for every pair of nodes in parallel. The first loop only runs $n - 2$ (instead of $n - 1$ times as I stated above) because you started with a matrix with all shortest path of length 1 already calculated (your base case is $t = 1$ instead of $t = 0$).

First loop of Extend-shortest-path is selecting the source node and following two loops are iterating through all edges (it assume the graph is complete) and relax the path that starts at $i$ ends in $j$ and use the edge $j-k$ next.



Old answer below (I miss-read the question and assumed it was asking  about Floyd-Warshall)

The key element that makes Floyd-Warshall (All-pairs-shortest-path) works is very different from Bellman-Fords key element.

Floyd-Warshall solves the following problem using dynamic programming. Calculate shortest path from every pair of vertex that use only first k vertex as intermediate hops. At the beginning it is the initial weight of existing edges and $\infty$ elsewhere (That is the base case).

Assume you have calculated such matrix for the first $k$ vertex and let's try to compute it for the first $k+1$ vertex. We should take into account paths that uses first $k$ vertex and goes through vertex $k+1$. Denote by $dp_k(u, v)$ the shortest path from $u$ to $v$ using first $k$ vertex, then it follows that the optimal path from $u$ to $v$ using $k+1$ vertex don't use this new vertex or indeed use it, so we get minimum from both values:

$$dp_{k+1}(u, v) = min(dp_k(u, v), dp(u, k) + dp(k, v))$$

In practice, the same table is reused on every step, so the algorithm only requires $O(n^2)$ spatial complexity. Temporal complexity is $O(n^3)$.