Besides SCAM in Matlab, there is also a slick online symbolic circuit analysis tool at CircuitNAV, which uses netlist files (from LTspice, Micro-Cap, TINA-TI, PSpice, etc) as the input and generates the algebraic solution for each circuit parameter.
CircuitNAV also provides a demo and a tutorial.
Doing that analytically will be really hard, even if you find an expression for \$ Q \$ there will be many local minima and you will have to test each of them to find the actual global one. But numerically you can find an approximate (or at least a local minimum). Do it like this using Octave (should be easy to adapt to Matlab),
pkg load control;
There is, you can solve it with the symbolic package. In Gnu Octave, I did as follows:
pkg load symbolic
syms s i1 i2 i3 i4
%Since the angles are in degrees, you'll need to use "cosd" and "sind" instead
eq_sys = [
(6*s*cosd(13)+12*sind(13))/(s^2+4) == (4000/(3*s)+2)*i1-4000*i2/(3*s)-2*i3;
1e-3 == (-4000/(3*s)+0.005)*i1+(4000/(3*s)+1e-3*s)*i2-1e-3*s*i4;...
Everyone has their own preferred simulation tool depending on what their background is and what they are simulating most often. You can get Almost all simulation tools to simulate any circuit(within reason) some just are a little easier in some scenarios. I would suggest you look at what simulation software the people around you or the people that you are ...