R=10;

z=real(zeros(201,201));

m=0;

n=0;

step = 0.1;

for x=-R:step:R

m = m + 1;

%x

for y=-sqrt(R*R - x*x):step:sqrt(R*R - x*x)

%y

n = int32(y / step) + R / step + 1;

%n

z(n, m) = real(sqrt(R*R - x*x - y*y));

endfor

endfor

%z

mesh(z);

另一种方法(from octave)：

function [xx, yy, zz] = sphere (varargin)

[hax, varargin, nargin] = __plt_get_axis_arg__ ("sphere", varargin{:});

if (nargin > 1)

print_usage ();

elseif (nargin == 1)

n = varargin{1};

else

n = 20;

endif

theta = linspace (0, 2*pi, n+1);

phi = linspace (-pi/2, pi/2, n+1);

[theta,phi] = meshgrid (theta, phi);

x = cos (phi) .* cos (theta);

y = cos (phi) .* sin (theta);

z = sin (phi);

if (nargout > 0)

xx = x;

yy = y;

zz = z;

else

oldfig = [];

if (! isempty (hax))

oldfig = get (0, "currentfigure");

endif

unwind_protect

hax = newplot (hax);

surf (x, y, z);

unwind_protect_cleanup

if (! isempty (oldfig))

set (0, "currentfigure", oldfig);

endif

end_unwind_protect

endif

endfunction