Your use of fitlm() is correct. Alternatively, if you don't need all the extra information provided by fitlm() and speed is a concern then you can use MATLAB mldivide (\.) to solve it more efficiently.
data = readtable('seyreg.csv');
X = [data.sst data.at data.vlm ones(size(data.sst))];
Y = data.sl;
coefficient = X\Y;
the coefficient are in same order as the column of X.