Get bmesh and list of vertices
Set  the points to a collection of coords.

Needed to change to the points from a collection of bmesh.types.BMVert to a collection of 3d coordinates.

points = [v.co.copy() for v in bm_points.verts]


May still need some tests, here is the result 



added some code such that if the grid is in edit mode, select the result.

import bpy
from mathutils import *
from mathutils.bvhtree import BVHTree
import bmesh


def are_inside(points, bm):
    """
    input: 
        points
        - a list of vectors (can also be tuples/lists)
        bm
        - a manifold bmesh with verts and (edge/faces) for which the 
          normals are calculated already. (add bm.normal_update() otherwise)
    returns:
        a list
        - a mask lists with True if the point is inside the bmesh, False otherwise
    """

    rpoints = []
    addp = rpoints.append
    bvh = BVHTree.FromBMesh(bm, epsilon=0.0001)

    # return points on polygons
    for point in points:
        fco, normal, _, _ = bvh.find_nearest(point)
        p2 = fco - Vector(point)
        v = p2.dot(normal)
        addp(not v &lt; 0.0)  # addp(v &gt;= 0.0) ?

    return rpoints



def test(grid, ob):
    print("test")
    mesh_points = grid.to_mesh()
    bm_points = bmesh.new()
    bm_points.from_mesh(mesh_points)
    points = [v.co.copy() for v in bm_points.verts]


    mesh = ob.data
    bm = bmesh.new()

    bm.from_mesh(ob.to_mesh())

    inside = are_inside(points, bm)
    if grid.mode == 'EDIT':
        # view the result
        bm = bmesh.from_edit_mesh(grid.data)

        for v, select in zip(bm.verts, inside):
            v.select = select
        bmesh.update_edit_mesh(grid.data)

scene = bpy.context.scene
test(scene.objects["Grid"], scene.objects["Suzanne"])


Consider using bmesh_from_object() or applying the global transform to the bmesh.  If using an edit mode bmesh and not wanting to commit changes bmesh.from_edit_mesh(me).copy() provides an unbound copy.

If only the transform needed then transforming a mesh copy with

me.transform(ob.matrix_world)


or using the operator, will provide global coordinates.

bmesh.ops.transform(bm,
        matrix=ob.matrix_world,
        verts=bm.verts)


Alternatively could put one in the object space of the other to make their coordinates equivalent.

The mesh provided from ob.to_mesh() is "not what it used to be"

How do I get a mesh data-block with modifiers and shape keys applied in Blender 2.8?

May be worth noting too that Suzanne has 3 separate parts.