def METHOD_NAME(fn, div_x=16, div_y=16, size_x=2, size_y=2):
    x, y = np.meshgrid(np.linspace(-size_x / 2, size_x / 2, div_x + 1),
                       np.linspace(-size_y / 2, size_y / 2, div_y + 1))
    z = fn(x, y)
    vertices = np.stack([x.flatten(), y.flatten(), z.flatten()]).T
    faces = np.array([[0, div_y + 1, div_y + 2, 1]]) + np.expand_dims(
        (np.expand_dims(np.arange(div_y), 0) + np.expand_dims(np.arange(div_x) * (div_y + 1), 1)).flatten(), -1)
    mesh = bpy.data.meshes.new('z_function_surface')
    mesh.from_pydata(vertices, [], faces)
    mesh.update()
    return mesh