arrayfun( @(A, B) roots([B^2+A^2+1, -12*B-10*A-8, 61], a, b, 'Uniform', 0)
Alternate approach that reaches the answer you want but with a different sequence of steps:
syms A B X
x = solve((X*A - 5)^2 + (X*B - 6)^2 + (X - 4)^2 - 16, X);
subs(x, {A, B}, {a, b})
This finds the generalized solution as a single equation, and then puts the actual values in, which is the reverse of the steps you asked.