Bones with Python: how to define kinematics from 4by4 homogeneous matrices?

Hi all,
I am a PhD student in robotics, and I am trying to use Blender to visualize the movement of a hand exoskeleton that I am working on. I ran into some difficulties trying to build the forward kinematics in Blender, and I am hoping that someone in this forum will be able to help me.

I have the CAD model of the exoskeleton, and I imported the single links as .stl files in Blender, each saved in their local coordinate system that I defined in my CAD program. Also from the CAD program I know the pose (orientation and translation) of each link respect to the one coming previously in the kinematic chain, defined as the relative transformation between local reference systems. Using this information I was able to write a Python script that assembles the exoskeleton in Blender in its rest position; however after I began trying to build an armature I started running intro problems.

The issue I am having is that my kinematic chain is defined as a sequence of rigid transformations that define the position of consecutive joints; however, as far as I was able to understand in Blender bones are defined by assigning the position of their head and tail and their roll angle, with the y axis assigned along the bone lenght. Basically what I need is a way to redefine the coordinate system for bone tails, which apparently is not a straightforward thing to do.

I spent a couple of days searching for a way to do this but I was not able to find a solution; the one that got closer to what I need is this post

However this is a hotfix for a single bone, and does not use python. I would like to be able to do this procedurally. Any suggestions on how to do this? If assigning the transformation is not possible it would be helpful to be able to obtain the homogeneous transformation of the bone tail respect to the previous coordinate system (that way I should be able to compensate for the difference with an offset).