Help setting up multivariate nonlinear equation problem

This might be a lot but I'll try to make it as concise and informative as possible.

Some background, part of a project I'm working on is to simulate joint motion (knee for example) where I know the initial and final body orientations, the shin and ankle in this example, then I also know the angles. So, the unknowns I have are the components of knee axes (flexion/extension [FE], abduction/adduction [AA], internal rotation [IR]), so 3 (x,y,z) vectors, but with additional equations/relationships I can define the AA and IR axes where with my coordinate system the AA axis is the cross product of FE and IR. So I'm left with just the FE axis as the unknown vector.

For my analysis, to get from initial to final body orientations I rotate the initial body segments about the axes using the known angles to end up at the final body orientation. I'm trying to find the 3 FE axis components that give the least error between the simulated final body orientations and the known ones for 100 time points. To get from initial to final involves at least 6 equations that involve rotating a vector (body segment) about an axis (joint axis) by an angle. Each equation depends on the one prior to it, whether it defines body segments or a new joint axis to use.

I did all the previous analysis using fsolve and lsqnonlin which was fine. I found that when I performed a similar analysis where I used the axis I previously solved with the known body segment orientations to solve for 3 angles about each axis at each time point I wouldn't get the angles I used in the previous analysis. So, now I'm trying to put another fsolve or lsqnonlin inside my previous one to ensure that the axis I solve for minimizes the error between the set of angles that correspond to the solved axis and the set of known angles because they aren't already equal for some reason. Another option is to make the angles unknown in my outside fsolve or lsqnonlin then set a condition to minimize the difference between the angles being solved for and what they should be. So I might be adding a nonlinear constraint into this as well.

Putting my angle solver inside my axis solver seems to work decently but it's very computationally expensive as it has taken approximately the past 10 hours to solve running my processor at 2.8 GHz. So, my main question here is 1) if there's a more efficient way to go about this that would give me accurate answers in a shorter time and 2) if there's a good way I can set this up to utilize parallel processing with my two solvers. I tried venturing into the parallel toolbox and other optimization functions but if it takes potentially up to 10 hours to figure out if I'm right or wrong I'd rather come here and try to get some advice.

I understand this is potentially a loaded and specific question but I'm more than willing to attach my code if requested. Thank you in advance.