Something like this should work:
x = 1:10;
y1 = rand(1, 10);
y2 = rand(1, 10);
DM = [x', ones(length(x),1)]; % Design Matrix
B1 = DM \ y1'; % Estimate Parameters Of ‘y1’ Fit
B2 = DM \ y2'; % Estimate Parameters Of ‘y2’ Fit
y1_fit = DM*B1; % Calculate Fit For ‘y1’
y2_fit = DM*B2; % Calculate Fit For ‘y2’
figure(1)
plot(x, y1, '.b', x, y2, '.r') % Plot Data
hold on
plot(x, y1_fit, '--b', x, y2_fit, '--r') % Plot Fitted Regression Lines
hold off
grid
