import pylab
from mpl_toolkits.mplot3d import Axes3D
import scipy.integrate as integrate
import numpy as np
from numpy import exp, sqrt, cos
a = 1.4
H = 10
def makeData():
    x = np.arange(0, 200, 2)
    y = np.arange(0, 100, 2)
    xgrid, ygrid = np.meshgrid(x, y)
    zgrid = integrate.quad(lambda k:2*(sqrt(k**2+a**2)*exp(-H*(k**2+a**2))*cos(xgrid*a*sqrt(k**2+a**2))*cos(ygrid*k*sqrt(k**2+a**2))),0, np.inf)
    return xgrid, ygrid, zgrid
x, y, z = makeData()
fig = pylab.figure()
axes = Axes3D(fig)
axes.plot_surface(x, y, z)
pylab.show()

 TypeError: only size-1 arrays can be converted to Python scalars

 #!/usr/bin/python3
import pylab
from mpl_toolkits.mplot3d import Axes3D
import scipy.integrate as integrate
import numpy as np
from numpy import exp, sqrt, cos
a = 1.4
H = 10
def f(x,y):
    return integrate.quad(lambda k:2*(sqrt(k**2+a**2)*exp(-H*(k**2+a**2))*cos(x*a*sqrt(k**2+a**2))*cos(y*k*sqrt(k**2+a**2))),0, np.inf)
def makeData():
    x = np.arange(0, 200, 2)
    y = np.arange(0, 100, 2)
    xgrid, ygrid = np.meshgrid(x, y)
    zgrid = np.array([f(x,y)[0] for x,y in zip(np.ravel(xgrid), np.ravel(ygrid))]).reshape(xgrid.shape)
    return xgrid, ygrid, zgrid
x, y, z = makeData()
fig = pylab.figure()
axes = Axes3D(fig)
axes.plot_surface(x, y, z)
pylab.show()

 <Figure size 640x480 with 1 Axes>

 %matplotlib inline

 zgrid = np.array([f(x,y)[0] for x,y in zip(np.ravel(xgrid), np.ravel(ygrid))]).reshape(xgrid.shape)

 #!/usr/bin/python3
import pylab
from mpl_toolkits.mplot3d import Axes3D
import scipy.integrate as integrate
import numpy as np
from numpy import exp, sqrt, cos
a = 1.4
H = 10
def f(x,y):
    return integrate.quad(lambda k:2*(sqrt(k**2+a**2)*exp(-H*(k**2+a**2))*cos(x*a*sqrt(k**2+a**2))*cos(y*k*sqrt(k**2+a**2))),0, np.inf)[0]
def makeData():
    x = np.arange(0, 200, 2)
    y = np.arange(0, 100, 2)
    xgrid, ygrid = np.meshgrid(x, y)
    zgrid = np.array([f(x,y) for x,y in zip(np.ravel(xgrid), np.ravel(ygrid))]).reshape(xgrid.shape)
    return xgrid, ygrid, zgrid
x, y, z = makeData()
fig = pylab.figure()
axes = Axes3D(fig)
axes.plot_surface(x, y, z)
pylab.show()

[0] в конце return integrate.quad(...

 >>> integrate.quad(lambda x:x,0,4)
(8.0, 8.881784197001252e-14)

 x = np.arange(0, 200, 2)

 >>> x = np.arange(0, 10, 4)
>>> print(x)
[0 4 8]
>>> y = np.arange(1, 10, 5)
>>> print(y)
[1 6]

 xgrid, ygrid = np.meshgrid(x, y)

 >>> print(xgrid)
[[0 4 8]
 [0 4 8]]
>>> print(ygrid)
[[1 1 1]
 [6 6 6]]

 zgrid = np.array([f(x,y) for x,y in zip(np.ravel(xgrid), np.ravel(ygrid))]).reshape(xgrid.shape)

 [f(x,y) for x,y in zip(np.ravel(xgrid), np.ravel(ygrid))]

 for x,y in zip(array1, array2)

 [f(x,y) for x,y in ..]

 np.array(список).reshape(xgrid.shape)

 ... .reshape(xgrid.shape)

 np.array(список).reshape((2,2))

 import scipy.integrate as integrate
import numpy as np
from numpy import exp, sqrt, cos
from mpmath import*
a = 1.4
H = 10
def opts0():
    return {'limit':1000}
def f(x,y):
    return integrate.nquad(lambda m,k:0.5*(k*exp(-H*k)*cos(a*(m+x)*sqrt(k))*besselj(0,k*sqrt(m**2+y**2))),[[0, np.inf],[-x,np.inf]],opts=[opts0, opts0, {}, {}],full_output=True)[0]

 f(1,1)

 TypeError: must be real number, not mpc

 C:\Users\all\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py:385: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  warnings.warn(msg, IntegrationWarning)
C:\Users\all\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py:385: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.
  warnings.warn(msg, IntegrationWarning)
C:\Users\all\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py:385: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)