How to find the concentrations of the acid-base equilibrium between sulfuric acid and sulfur dioxide for a known pH using Octave?
Your code seems to be calculating the differences between a guessed concentration vector x and a calculated concentration vector.  These differences in concentrations between guessed and calculated are y.  For example, the guessed pH is apparently x(1), and the calculated pH is 10^(-pH).  The difference between them is y(1).

Thus, its no surprise that the y values are sometimes negative.  To solve the problem you need to feed your function to a numerical solver such as fsolve.  That function will take your guessed solution x, other parameters such as pH, the various pKas, and M, and your function phCalc as input, and return the value for x which makes all the y variables zero.  That solution is the one you want.

On the chemistry side, you mention sulfur dioxide, but I don't see any equations for the solubility of sulfur dioxide gas in water.  Additionally, sulfur dioxide is reduced relative to sulfuric acid, so there would need to be an additional reaction, and parameter (extent of reaction or degree of reduction or something) that describes the relative amounts of reduced sulfur (as sulfur dioxide, bisulfite, or sulfite) vs. oxidized sulfur (as sulfate, bisulfate, sulfuric acid, or sulfur trioxide).

There is also the issue of charge balance: what concentration of inert cations (e.g. sodium or potassium or whatever) are you assuming?  Charge balance must hold in bulk aqueous solution.  If you don't have any sodium / potassium etc., then you should include the self ionization of water as an equation, hydroxide as a species, and y(7) = x(1) - x(2) - 2*x(3) - x(4) - x(OH) as an equation.  If you do have sodium/potassium/whatever at significant concentration, you can neglect hydroxide and just go with the equation y(7) = x(1) - x(2) - 2*x(3) - x(4) + x(Na) where x(Na) is the cation concentration, which I have assumed is sodium.

#