Deriving a formula for $b_k$ is straightforward. You can multiply each side by $(x+k)^2$, simplify, and set $x = -k$ to see that$$
b_k = \left(\frac{1}{(-k)(-(k-1))\cdots (-1)(1)(2)\cdots(n-k)}\right)^2 = \frac{1}{\left(k!(n-k)!\right)^2},
$$
in agreement with your calculations.
Deriving a formula for $a_k$ is more challenging. The rising and falling factorial functions may be expanded as polynomials whose coefficients are the unsigned Stirling numbers of the first kind:$$
x^{(n)} = x(x+1)\cdots (x+n-1) = \sum_{j=0}^{n}{n \brack j} x^j,\;\text{and}\\
(x)_n = x(x-1)\cdots (x-(n-1)) = \sum_{j=0}^{n}(-1)^{n-j}{n \brack j} x^j.
$$
The goal is to find coefficients of an expansion of $\left(x^{(n+1)}\right)^{-2}$ in powers of $y \doteq x+k$. Therefore we write
$$
x^{(n+1)}= (y-k)(y-(k-1)) \cdots (y-1) y (y+1) \cdots (y+n-k) = \frac{y^{(n-k+1)}(y)_{k+1}}{y}.
$$
Expanding this in powers of $y$ yields
$$
x^{(n+1)} = \frac{1}{y}\left({n-k+1 \brack 0} + {n-k+1 \brack 1} y+ {n-k+1 \brack 2} y^2 + O(y^3)\right) \times
(-1)^k \left(-{k+1 \brack 0} + {k+1 \brack 1} y - {k+1 \brack 2} y^2 + O(y^3)\right)
$$
The formulas for the Stirling numbers involved are, for all $n \ge 0$, $$
{n+1 \brack 0} = 0, \quad
{n+1 \brack 1} = n!, \quad \text{and} \quad
{n+1 \brack 2} = n! H_n,
$$
where $H_n$ is the $n^\mathrm{th}$ harmonic number. Therefore,
$$
x^{(n+1)} = (-1)^k k!(n-k)! \, y \left(1 + (H_{n-k} - H_k)y + O(y^2)\right).
$$
Raising this to the power $-2$ yields
$$
\left(x^{(n+1)}\right)^{-2} = \frac{1}{\left(k!(n-k)! \right)^2}\, \left(y^{-2} + 2 (H_k - H_{n-k}) y^{-1} + O(1)\right).
$$
From this we may read off the formula for $b_k$ given above and this formula for $a_k$: $$
a_k = 2\frac{H_k - H_{n-k}}{\left(k!(n-k)! \right)^2}.
$$
Here's an example for $n = 3$. The values of $b_k$ are $$b_0 = b_3 = \frac{1}{(0!3!)^2} = \frac{1}{36} \;\text{and}\; b_1 = b_2 = \frac{1}{(1!2!)^2} = \frac{1}{4}.
$$
The values of $a_k$ are $2(H_k - H_{n-k})b_k$.
The first few harmonic numbers are$$
H_0 = 0,\;H_1 = \frac{1}{1} = 1,\;H_2 = \frac{1}{1} + \frac{1}{2} = \frac{3}{2},\;\text{and}\;H_3 = \frac{1}{1} + \frac{1}{2} + \frac{1}{3} = \frac{11}{6}.
$$
Therefore $$
a_3 = 2(H_3 - H_0)b_3 = 2\left(\frac{11}{6} - 0\right)b_3 = \frac{11}{3}\cdot\frac{1}{36} = \frac{11}{108},\\
a_2 = 2(H_2 - H_1)b_2 = 2\left(\frac{3}{2} - 1\right)b_2 = \frac{1}{1}\cdot\frac{1}{4} = \frac{1}{4},
$$
whereas $a_1 = -a_2$ and $a_0 = -a_3$ because of the (anti)symmetry in the formula for $a_k$. Therefore $$
\frac{1}{\left(x(x+1)(x+2)(x+3)\right)^2} = -\frac{11}{108} \frac{1}{x} - \frac{1}{4} \frac{1}{x+1} + \frac{1}{4} \frac{1}{x+2} + \frac{11}{108} \frac{1}{x+3} + \frac{1}{36} \frac{1}{x^2} + \frac{1}{4} \frac{1}{(x+1)^2}+ \frac{1}{4} \frac{1}{(x+2)^2}+ \frac{1}{36} \frac{1}{(x+3)^2}.
$$