ここでは、

    b0                   = p0
    b0 +  b1 + ... +   bl = p1
    b0 + 2b1 + ... + 2lbl = p2
               ...
    b0 + lb1 + ... + llbl = pl
               ...
    b0 + nb1 + ... + nlbl = pn
 

を辺々引いて b_i が解きやすくなるように右三角形にします。 各行から一つ上の行を引く、という操作を i 行目なら i - 1回行います。

引き算に b_i は関係無いので、係数だけ考えればよいです。 各行の係数を一括して考えるために各式を次のようにします。

    1 + ix/1! + i2x2/2! + ... = yi         (1)
 

p_i も係数が1の変数と考えます。
例えば、元の3番目の式から2番目の式を引くと、

    b1 + 3b2 + ... + (2l - 1)bl = p3 - p2
 

となりますが、(1)の引き算も係数だけ見れば同じになります。

    x/1! + 3x2/2! + ... + (2l - 1)xl/l! + ... = y3 - y2
 

項を後ろに追加して無限級数にしていますが、影響はありません。 こうするのは、(1)は次のように簡単になります。

    eix = yi
 

ここで、

    fi(x) = eix
    fi0(x) = fi(x)
    fij+1(x) = fij(x) - fi-1j(x)
 

とおきます。計算していくと、

    fi1(x) = (ex - 1)e(i-1)x
    fi2(x) = (ex - 1)2e(i-2)x
            …
    fii(x) = (ex - 1)i
 

となることが帰納的にわかります。

    fii(x) = sum(iCkekx(-1)i-k, k=0,…,i)
    f(j)ii(x) = sum(iCkekx(-1)i-kkj, k=0,…,i)
 

結局 b_j の各係数 a_ij は、

    aij = f(j)ii(0) = sum(iCk(-1)i-kkj, k=0,…,i)
 

また、

    fii(x) = (ex - 1)i
           = (x + x2/2! + …)i
           = xi + …
 

だから、

    aij = 0  (j < i)
        = i! (j = i)
 

右辺も同様に求まります。
もっとも、本当にこれに基づいて計算するほうが、 一つ一つ引き算を計算するより速いかどうかは確かめていません。