4841.sagews

R.<p0,p1,beta,gamma> = PolynomialRing(ZZ,4)

P=matrix(R,2,[1-p0,p0,1-p1,p1])

B=matrix(R,2,[beta,1,1,gamma])

M=P*B*P.transpose()

(M[0,0]-M[0,1]) == -(p0-p1)*((beta-1)*(1-p0)+p0*(1-gamma))

True