For https://blender.stackexchange.com/a/279753/88681

material_switch_anim.py

import bpy

#### SETUP MESH AND MATERIALS

# Create mesh

bpy.ops.mesh.primitive_cube_add()
ob = bpy.context.object

# Create materials

def make_material(name, color):
    mat = bpy.data.materials.new(name)
    mat.use_nodes = True
    node = mat.node_tree.nodes["Principled BSDF"]
    node.inputs["Base Color"].default_value = color + (1,)
    return mat

red_mat = make_material("Red", (1, 0, 0))
blue_mat = make_material("Blue", (0, 0, 1))

# Add materials to slots

ob.data.materials.append(red_mat)    # Red = slot 0
ob.data.materials.append(blue_mat)   # Blue = slot 1


#### GEOMETRY NODES

# Create Geometry Nodes input/output

group = bpy.data.node_groups.new("Geometry Nodes", 'GeometryNodeTree')
group.inputs.new('NodeSocketGeometry', "Geometry")
group.outputs.new('NodeSocketGeometry', "Geometry")
input_node = group.nodes.new('NodeGroupInput')
output_node = group.nodes.new('NodeGroupOutput')
output_node.is_active_output = True

# Add Set Material Index node

set_node = group.nodes.new('GeometryNodeSetMaterialIndex')
group.links.new(input_node.outputs[0], set_node.inputs[0])
group.links.new(set_node.outputs[0], output_node.inputs[0])

# Position nodes for aesthetics

input_node.location = -160, 40
set_node.location = 12, 68
output_node.location = 190, 50

# Add Geometry Node modifier

mod = ob.modifiers.new("Material Switcher", 'NODES')
mod.node_group = group


#### ANIMATION

anim = [
    # frame,  material slot
    0,        0,
    5,        1,
    10,       0,
    11,       1,
    12,       0,
    15,       1,
]

# Insert keyframes with keyframe_insert

for i in range(0, len(anim), 2):
    frame, slot = anim[i], anim[i+1]

    set_node.inputs[2].default_value = slot
    group.keyframe_insert(
        'nodes["Set Material Index"].inputs[2].default_value',
        frame=frame,
    )

# Turn off interpolation

action = group.animation_data.action
kps = action.fcurves[0].keyframe_points
for kp in kps:
    kp.interpolation = 'CONSTANT'