def METHOD_NAME(tape, wire_order=None):
    """Defines how matrix works if applied to a tape containing multiple operations."""
    if not tape.wires:
        raise qml.operation.MatrixUndefinedError
    params = tape.get_parameters(trainable_only=False)
    interface = qml.math.get_interface(*params)
    wire_order = wire_order or tape.wires
    # initialize the unitary matrix
    if len(tape.operations) == 0:
        result = qml.math.eye(2 ** len(wire_order), like=interface)
    else:
        result = matrix(tape.operations[0], wire_order=wire_order)
    for op in tape.operations[1:]:
        U = matrix(op, wire_order=wire_order)
        # Coerce the matrices U and result and use matrix multiplication. Broadcasted axes
        # are handled correctly automatically by ``matmul`` (See e.g. NumPy documentation)
        result = qml.math.matmul(*qml.math.coerce([U, result], like=interface), like=interface)
    return result