close all;
clear all;

x0 = [0 0 10];
x1 = [1 0 0];
x2 = [0 1 0];
x3 = [0 0 0.5];
x4 = [0 0 -0.25];

N = 10;
M = 2*N+1;

% define a target arrays
p = zeros( M,M,3 );

ri = -N:N;
sj = -N:N;
% compute the position at each (r,s) pair
for i =1:M
  for j = 1:M
    r = ri(i);
    s = sj(j);
    p(i,j,:) = x0 + r*x1 +  r^2*x3 + s*x2 + s^2*x4;
  end
end

% make a visualisation
figure;
surf( p(:,:,1), p(:,:,2), p(:,:,3) );


figure;
hold on
surf( p(:,:,1), p(:,:,2), p(:,:,3) );
surf( p(:,:,1), p(:,:,2), 0*p(:,:,3) );
view(3);