If there were no interaction between x and y (c3=0), then c1 and c2 could be easily found by the backslash operator in MATLAB in a least square way (A\b)
Backslash should still work for least squares estimation,
x=x(:); y=y(:); z=z(:);
c = [x,y,x.*y].\z; %c=[c1;c2;c3]
although, if you have measurement noise in x and y (regardless of c3=0), some would say you should be using total least squares instead of ordinary least squares,
x=x(:); y=y(:); z=z(:);
A=[x,y,x.*y,-z];
A=bsxfun(@minus,A,mean(A,1));
[~,~,V]=svd(A,0);
c=V(1:end-1,end)./V(end);
