How to plot a perfectly smooth sphere

I am trying to plot a perfectly smooth sphere in python using matplotlib. I have been using the following code:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

fig = plt.figure(1)    
ax = fig.add_subplot(111, projection='3d')

u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))

ax.plot_surface(x, y, z, linewidth=0.0)

plt.show()

The plot I obtained is attached below:

The figure generated consists of patches of rectangular surfaces on it. Is it possible to smoothen out the boundaries or make it indistinguishable and to make a perfectly smooth sphere?

Matplotlib plots 3d surfaces by breaking them down into small sub-polygons of equal colour, as is explained in the documentation, hence your result is not really a surprise. In order to get a smoother surface, you have to provide more data points. There is, however, a small twist, which is that plot_surface() may not use all the data you provide. This is controlled with the cstride and rstride keywords. How the defaults are calculated is not quite clear to me, but below a little example that demonstrates the effect:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

fig,axes = plt.subplots(ncols=2,nrows=2,subplot_kw=dict(projection='3d'))

N=50
stride=2
ax = axes[0,0]
u = np.linspace(0, 2 * np.pi, N)
v = np.linspace(0, np.pi, N)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, linewidth=0.0, cstride=stride, rstride=stride)
ax.set_title('{0}x{0} data points, stride={1}'.format(N,stride))

N=50
stride=1
ax = axes[0,1]
u = np.linspace(0, 2 * np.pi, N)
v = np.linspace(0, np.pi, N)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, linewidth=0.0, cstride=stride, rstride=stride)
ax.set_title('{0}x{0} data points, stride={1}'.format(N,stride))

N=200
stride=2
ax = axes[1,0]
u = np.linspace(0, 2 * np.pi, N)
v = np.linspace(0, np.pi, N)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, linewidth=0.0, cstride=stride, rstride=stride)
ax.set_title('{0}x{0} data points, stride={1}'.format(N,stride))

N=200
stride=1
ax = axes[1,1]
u = np.linspace(0, 2 * np.pi, N)
v = np.linspace(0, np.pi, N)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, linewidth=0.0, cstride=stride, rstride=stride)
ax.set_title('{0}x{0} data points, stride={1}'.format(N,stride))

plt.show()

The resulting figure looks like this:

As you can see, the outcome of the plot is sensitive to both the density of your data and the stride keywords. Be careful though with the amount of data you provide -- plot_surface() can take up a considerable amount of time to provide a result. Hope this helps.