function pup()
    [x y] = meshgrid(-1 : 0.02 : 1, -1 : 0.02 : 1);
    d = x.^2 + y.^2;
    i1 = d < 0.25;
    i2 = ~i1 & d < 1;
    z = i1 * sqrt(0.75) + i2 .* sqrt(1 - d);
    surf(x, y, z);
end