Per the doc, the modelspec for a quadratic is 'purequadratic' or, with only the one variable, 'quadratic' will end up with the same model as there are no variables to do the interaction terms between.
You'll have to call fitlm once for each set of X,Y; the function isn't vectorized to compute multiple responses in one call. You can do this most effeciently coding-wise by having the data in an array and address the columns programmatically in a for...end loop rather than writing the code five-times over.
