 Solution

Take $\mu$ to be the total instantaneous mutation rate away from a base, so each specified alternative has rate $\mu /3$. Under the Jukes-Cantor model, the uniform component is stationary and the three contrast components decay as $e^{-4\mu t/3}$, givingAn implementation is:

def jukes(a, b, t, mu):
    decay = np.exp(-4*mu*t/3)
    return .25 + .75*decay if a == b else .25*(1-decay)

print(f"{jukes(2,2,1,.5):.3f} {jukes(2,3,1,.5):.3f}")

The two probabilities are $0.635$ and $0.122$. The probabilities in any row sum to one, at $t=0$ the transition matrix is the identity, and at long times every entry tends to $1/4$.

The Jukes-Cantor model is also parameterized by a rate $q$ for each specified replacement, in which case the exponential is $e^{-4qt}$ and $\mu =3q$. The total-rate convention is used consistently below; it yields an optimum within the stated search interval, whereas the corresponding per-replacement rate is about $0.0522$.