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

remove_originals = True
remove_doubles = True

context = bpy.context
scene = context.scene
# by selection
cubes = context.selected_objects
# by name
#cubes = [o for o in scene.objects if o.name.startswith("Cube")]

bm = bmesh.new()
for o in cubes:
    bm.from_mesh(o.data)
    bmesh.ops.transform(bm,
            verts=bm.verts[-8:],
            matrix=o.matrix_world
            )
if remove_originals:
    bpy.ops.object.delete()

bvhtree = BVHTree().FromBMesh(bm, epsilon=1e-7)
faces = bm.faces[:]

remove = list()
while faces:
    f = faces.pop()
    pair = bvhtree.find_nearest_range(f.calc_center_median(), 1e-4)
    if len(pair) > 2:
        # mark face for removal
        remove.extend(p[2] for p in pair)

bm.faces.ensure_lookup_table()
bmesh.ops.delete(bm,
        geom=[bm.faces[i] for i in set(remove)],
        context='FACES_KEEP_BOUNDARY',
        )
if remove_doubles:
    bmesh.ops.remove_doubles(bm,
            verts=bm.verts,
            dist=1e-7)

me = bpy.data.meshes.new("CubeBoundaryFaces")
bm.to_mesh(me)
ob = bpy.data.objects.new("CubeBoundaryFaces", me)

scene.collection.objects.link(ob)
context.view_layer.objects.active = ob
ob.select_set(True)