If they really are multivariate polynomials (e.g., no sin() or gamma() or integrals, no sqrt() or x^n when n is not 0 or a positive integer), then the \ operator might be appropriate, especially if it is only the right-hand sides that are varying (in which case they can all be done at the same time.)
If the coefficients are changing each time, I would consider using symbolic algebra to solve() for the variables, and then turn the resulting formula into MATLAB code (matlabFunction() might help.) Any RootOf() symbolic polynomials can be converted in to roots() calls.
