Try this:
DM = [x(:).^2 ones(size(x(:)))]; % Design Matrix
B = DM \ y(:); % Parameters
yfit = DM * B; % Calculated Fit
SStot = sum((y-mean(y)).^2); % Total Sum-Of-Squares
SSres = sum((y(:)-yfit(:)).^2); % Residual Sum-Of-Squares
Rsq = 1-SSres/SStot; % R^2
To plot it:
figure
plot(x, y, 'p')
hold on
plot(x, yfit, '-r')
hold off
grid
