MATLAB coding problem with least squares surface fit

[EDIT: Thu May 19 17:56:20 UTC 2011 - Reformat - MKF]

Good (Proper Greeting of the Day All),

Again, I am very new to the MATLAB programming thing but I think I am getting better (yay!!) so any advice or help would be much appreciated. I am currently working on a program that will take data (x,y,z) and generate a surface from them using least squares approximation method. I am pretty sure this feature exists in one of MATLAB's amazing toolbars but I decided to take it upon myself and try to make my own. I want to plot the surface as well as give the error output. This is what I have so far but it says error:

??? Error using ==> mldivide
Matrix dimensions must agree.

Here is my coding:

function f=planefit(x,y,z)

%PLANEFIT(x, y, z) fits a least-squares polynomial through
%a set of data where x and y and z are the coordinates.

n = length(x); 
S1 = sum(x); S2 = sum(y); S3 = sum(z);
S4 = sum(x.*x); S5 = sum(x.*y); S6 = sum(y.*y);
S7 = sum(z.*x); S8 = sum(z.*y);

%Solving Matrix
A = [S4 S5 S1; S5 S6 S2; S1 S2 n];
B = [S7 S8 S3]; 
C = A\B;

%Surface Polynomial
for k=1:n
z = C(1)*x(k) + C(2)*y(k) + C(3); 
end

%Error Calculation
Error(k) = abs(z(k)-z);
Error = sum(Error.*Error);

%data visualization
[X Y}=meshgrid(x,y)
figure
surf(X,Y,z,'r')  %The Least Squares Surface
hold on
plot3(x, y,z 'o')  %The Original Data Points
grid on

EDIT After trying Sean de's suggestion.

Here is the new error I now have:

Error =

       0

??? Attempted to access z(2); index out of bounds because numel(z)=1.

Error in ==> planefit at 19 Error(k) = abs(z(k)-z)

Here is my current updated code:

function f=planefit(x,y,z)
%PLANEFIT(x, y, z) fits a least-squares polynomial through
%a set of data where x and y and z are the coordinates.
n = length(x);
S1 = sum(x); S2 = sum(y); S3 = sum(z);
S4 = sum(x.*x); S5 = sum(x.*y); S6 = sum(y.*y);
S7 = sum(z.*x); S8 = sum(z.*y);
%Solving Matrix
A = [S4 S5 S1; S5 S6 S2; S1 S2 n];
B = [S7; S8; S3];
C = A\B;
%Surface Polynomial
for k=1:n
    z = C(1)*x(k) + C(2)*y(k) + C(3);
    Error(k) = abs(z(k)-z)
end
%Error Calculation
Error = sum(Error.*Error)
%data visualization
[X Y]=meshgrid(x,y); figure
surf(X,Y,z,'r') %The Least Squares Surface
hold on
plot3(x, y,z,'o') %The Original Data Points
grid on
plot3(x, y,z,'o') %The Original Data Points
grid on

Answer 1

[Sean de Wolski]

[EDIT: Thu May 19 18:44:56 UTC 2011 - Merged Comments - MKF]

Should B be a column vector? Either:

B = [S7;S8;S3];

or

B = [S7 S8 S3].';

EDIT

for k = 1:n 
   z(k) = stuff 
%more stuff 
end

right now you're overwriting z on each loop iteration.

[Melissa]

I did not think B had to be but I will try it.

[Sean de Wolski]

Inner matrix dimensions must agree:

[3 3][1 3] does not work but [3 3][3 1] does.

[Melissa]

that error went away but now I have generated more...thank you!

Answer 2

[Tejaswini Gautham]

What is the matlab built-in function to do this?