I would do this numerically, using fsolve. It requires a slight re-write of your equations to make them all implicit.
Example —
syms Qa_1 Q1_1 Q2_1 Q3_1 Q4_1 real
eqn1 = (Qa_1 - (Q1_1 + Q2_1 + Q3_1 + Q4_1));
eqn2 = (Qa_1^2/60.51 + Q1_1^2/0.8616 - (1.035/Q1_1 + 24.3/Qa_1));
eqn3 = (Qa_1^2/60.51 + Q2_1^2/1.346 - (1.321/Q2_1 + 16.57/Qa_1));
eqn4 = (Qa_1^2/60.51 + Q3_1^2/1.346 - (1.236/Q3_1 + 8.873/Qa_1));
eqn5 = (Qa_1^2/60.51 + Q4_1^2/1.346 - (1.044/Q4_1 + 1.619/Qa_1));
Eqnsfcn = matlabFunction([eqn1, eqn2, eqn3, eqn4, eqn5], 'Vars',{[Qa_1, Q1_1, Q2_1, Q3_1, Q4_1]});
B0 = rand(1,5)*100;
[B,fval] = fsolve(Eqnsfcn, B0)
This was almost instantaneous. There are likely multiple roots, so experiment with different initial parameter estimates (here ‘B0’).
EDIT — This version makes it easier to track the individual variable names:
Eqnsfcn = matlabFunction([eqn1, eqn2, eqn3, eqn4, eqn5], 'Vars',{Qa_1, Q1_1, Q2_1, Q3_1, Q4_1});
B0 = rand(1,5)*100;
[B,fval] = fsolve(@(b)Eqnsfcn(b(1),b(2),b(3),b(4),b(5)), B0)
