How to combine multiple cubes with double faces?
BMesh script.



Wrote the script for another robloxic type question, buggered if I can find it maybe didn't post an answer. Also note script is for blender 2.8, whereas your question image is clearly 2.79. If you are just starting out scripting in blender, my advice is to start using 2.8

Any way, here is a script to run in object mode with all the cubes you want to join selected.

Loads all selected objects into a global bmesh.  Make sure only cubes are selected, haven't tested for it as script expects selected objects to only be 8 vertex cubes.  Also added commented out code to select cubes based on their name starting with "Cube"

Tests if two faces have same center using a bvh tree loaded with the bmesh. The find_nearest_range method finds each face twice, hence 4 matches being a pair, and 2 matches being no match

A further more robust test would be checking opposite facing normals, and to be sure to be sure the vert coordinates of the face themselves.

If remove_originals is set to True then the original cube objects are removed.

If remove_doubles is set to True it will remove doubles.

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) &gt; 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)


Note, would be for quicker,   as pointed out by subd guru @RobinBetts, to simply join, remove doubles then  remove interior faces, done via a script using the appropriate operators and object mode for them.

Another take would be to use the set of face indices to remove the faces from the original cubes without joining them. Since the compound bmesh faces will be in same order as the originals.

for i in set(remove):
    ob = cubes[i // 6] # 6 faces per cube
    f = ob.data.polygons[i % 6]