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

def f(x,y):
     r = np.sqrt(x**2 + y**2)
     result = np.sin(r)/r
     result[r==0] = 1.0
     return result

x_ticks = np.linspace(-10,10,51)
y_ticks = np.linspace(-10,10,51)

x,y = np.meshgrid(x_ticks,y_ticks)
z = f(x,y)
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(x,y,z,rstride = 1,cstride = 1,cmap = cm.Y1GnBu_r)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')