Using Python to Copy Rotation of Pose Bones Description Of What I Am Doing

I am building a rig that uses both IK and FK controls. I want this rig to have a script which will set the rotations of my FK control bones to match the pose of the IK bones (similar to rigify).

My method is to do the following:


Name the IK and FK bones similarly
Use bpy to collect the FK bones I desire and find their corresponding IK counterparts
Set the pose matrices of the FK bones to their corresponding IK bones


Here is the code I use:

# ... code will define variables for armature and id for IK-chain   

# move all FK bones to the current IK orientation
for bone in armature.bones:
    if "MCH-finger-fk." + IK_id in bone.name:
        ik_name = bone.name[:11] + "ik" + bone.name[13:]
        ik_bone = armature.bones[ik_name]
        bone.matrix = ik_bone.matrix


Results...

This almost works! The first bone in the FK sequence matches the corresponding IK bone perfectly. The remaining bones do not, but if I run it again, the next bone in the sequence will match! That is, the second time I run the script, the second bone in the FK chain will match its corresponding IK bone. 

The pictures below will hopefully explain what I mean.

Before Script Execution



Script Executed Once



Script Executed Twice

Alas, I cannot have more than two links, due to my low reputation. Imagine moving the first separated bone to where it looks like it should go (the remaining bones move relative) in the previous image.

Observation

I assume this bug in my script has something to do with the fact that only the second time I execute the script, the bone data takes the parent matrix into account (since the pose is now in place).

My Question

Is there a nicer way to achieve what I am trying? In essence, I want to copy the world rotations of each bone. I tried looking for the source code of the COPY_ROTATION constraint, but I was unsuccessful.
See this section of the api documentation. Setting PoseBone.matrix probably does the calculations required for PoseBone.matrix_basis but does not update the value of the underlying DNA world matrix property, which is wierd. Here is a script which uses the scene.update method: 

import bpy

class BoneMatrixCopier:
    def __init__(self, obj, scene):
        self.obj = obj
        self.scene = scene
        self.pose_bones = obj.pose.bones

    def update(self):
        self.obj.update_tag({'OBJECT'})
        self.scene.update()

    def get_ik_bone_name(self, name):
        return name.replace("Bone", "BoneIK")

    def copy_matrix(self, bone):
        self.update()
        bone.matrix = self.pose_bones[
            self.get_ik_bone_name(bone.name)
            ].matrix

    def copy_matrices(self, root_bone_name):
        root_bone = self.pose_bones[root_bone_name]

        self.copy_matrix(root_bone)
        for bone in root_bone.children_recursive:
            self.copy_matrix(bone)

bmc = BoneMatrixCopier(bpy.data.objects['Armature'], bpy.context.scene)
bmc.copy_matrices('Bone')