syms q1 q2 q3 q12 b1 b2 a1 a2 % eqn1 = q1 * a1 + (q3 + (q12/2)) *a2 == -2*q1*b1 - q12*b2 % % eqn2 = (q3 + (q12/2)) * a1 + q2 *a2 == -2*q2*b2 - q12*b1 eqn1 = q1 * a1 + (q3 + (q12/2)) *a2 == q1*b1 + (q12*b2/2) eqn2 = (q3 + (q12/2)) * a1 + q2 *a2 == q2*b2 + (q12*b1/2) [A,B] = equationsToMatrix([eqn1, eqn2], [a1,a2]) X = linsolve(A,B) % check special case fprintf('special case q12 = 0 ') subs(X,q12,0)