The method you used does not include an intercept, so your fit is forced to go through the origin. Try this instead (which is also on that documentation page).
y = [-38.34; -42.48; -52.71; -49.72];
x = [69.482; 60.645; 37.093; 42.889];
X = [ones(size(x)), x];
b1 = X\y;
yreg = X*b1;
scatter(x,y)
hold on
plot(x,yreg)
