from sympy import linsolve, symbols, Matrix, simplify

A, B, f, n, t, b, l, r = symbols('A B f n t b l r')

S = Matrix([
 [2/(r-l),       0,      0,   0],
 [0      , 2/(t-b),      0,   0],
 [0      ,       0,2/(f-n),   0],
 [0      ,       0,      0,   1]   
])

T = Matrix([
 [1, 0, 0,   -(r+l)/2],
 [0, 1, 0,   -(t+b)/2],
 [0, 0, 1,   -(-n-f)/2],
 [0, 0, 0,   1]   
])

P = Matrix([
 [n, 0, 0,   0],
 [0, n, 0,   0],
 [0, 0, A,   B],
 [0, 0, -1,  0]   
])

F = Matrix([
 [1, 0, 0,  0],
 [0, 1, 0,  0],
 [0, 0, -1, 0],
 [0, 0, 0,  1]   
])

O = simplify(S*T)

display(O)

$\displaystyle \left[\begin{matrix}\frac{2}{- l + r} & 0 & 0 & \frac{l + r}{l - r}\\0 & \frac{2}{- b + t} & 0 & \frac{b + t}{b - t}\\0 & 0 & \frac{2}{f - n} & \frac{f + n}{f - n}\\0 & 0 & 0 & 1\end{matrix}\right]$

\begin{equation} \begin{bmatrix} n & 0 & 0 & 0 \\ 0 & n & 0 & 0 \\ 0 & 0 & A & B \\ 0 & 0 & -1 & 0 \end{bmatrix} \begin{bmatrix} x \\ y \\ -f \\ 1 \\ \end{bmatrix} = \begin{bmatrix} nx \\ ny \\ -f \cdot f \\ f \\ \end{bmatrix} \sim \begin{bmatrix} \frac{nx}{f} \\ \frac{ny}{f} \\ -f \\ 1 \\ \end{bmatrix} \end{equation}\begin{equation} \begin{bmatrix} n & 0 & 0 & 0 \\ 0 & n & 0 & 0 \\ 0 & 0 & A & B \\ 0 & 0 & -1 & 0 \end{bmatrix} \begin{bmatrix} x \\ y \\ -n \\ 1 \\ \end{bmatrix} = \begin{bmatrix} nx \\ ny \\ -n \cdot n \\ n \\ \end{bmatrix} \sim \begin{bmatrix} x \\ y \\ -n \\ 1 \\ \end{bmatrix} \end{equation}

#eqns = [-A * n + B + n*f, -A * f + B + n*f] # flip frustrum
eqns = [-A * n + B + n*n, -A * f + B + f*f]  # no flip
linsolve(eqns, A, B)

$\displaystyle \left\{\left( f + n, \ f n\right)\right\}$

P = P.subs(A, f+n).subs(B, f*n)

# execute full pipeline
simplify(F*O*P)

$\displaystyle \left[\begin{matrix}\frac{2 n}{- l + r} & 0 & \frac{l + r}{- l + r} & 0\\0 & \frac{2 n}{- b + t} & \frac{b + t}{- b + t} & 0\\0 & 0 & - \frac{f + n}{f - n} & - \frac{2 f n}{f - n}\\0 & 0 & -1 & 0\end{matrix}\right]$