I’m attempting to model the motion of a figure skater, and currently have to solve this differential equation, numerically, for an angular velocity vector, (in the skater’s reference frame), and I’m having trouble understanding why this differential equation is failing to solve.

Here is the equation:

```
sol = NDSolve({(Li'(
t) + (ScriptCapitalI)'(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(
t)} + (ScriptCapitalI)(
t).{(Omega)x'(t), (Omega)y'(t), (Omega)z'(t)} +
Cross({(Omega)x(t), (Omega)y(t), (Omega)z(t)},
Li(t) + (ScriptCapitalI)(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(t)})).ê0 ==
0, (Li'(t) + (ScriptCapitalI)'(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(
t)} + (ScriptCapitalI)(
t).{(Omega)x'(t), (Omega)y'(t), (Omega)z'(t)} +
Cross({(Omega)x(t), (Omega)y(t), (Omega)z(t)},
Li(t) + (ScriptCapitalI)(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(t)})).ê1 ==
0, (Li'(t) + (ScriptCapitalI)'(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(
t)} + (ScriptCapitalI)(
t).{(Omega)x'(t), (Omega)y'(t), (Omega)z'(t)} +
Cross({(Omega)x(t), (Omega)y(t), (Omega)z(t)},
Li(t) + (ScriptCapitalI)(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(t)})).ê2 ==
0, (Omega)x(0) == (Omega)x0, (Omega)y(
0) == (Omega)y0, (Omega)z(0) == (Omega)z0}, {(Omega)x(
t), (Omega)y(t), (Omega)z(t)}, {t, 1, 2},
Method -> {"Automatic", {"EquationSimplification" -> "Residual"}});
```

From this full code:

```
Remove("Global`*");
(*describe instantaneos orientation of each piece, angles relative to
adjoining piece that is closer to center of body, when all angles are
set to zero the body is standing upright with arms hanging at sides*)
(Theta)1(t_) = 0; (Phi)1(t_) = 0; (Psi)1(t_) = 0;
(Theta)2(t_) = 0; (Phi)2(t_) = 0; (Psi)2(t_) = 0;
(Theta)3(t_) = 0; (Phi)3(t_) = 0; (Psi)3(t_) = 0;
(Theta)4(t_) = 0; (Phi)4(t_) = 0; (Psi)4(t_) = 0;
(Theta)5(t_) = 0; (Phi)5(t_) = 0; (Psi)5(t_) = 0;
(Theta)6(t_) = 0; (Phi)6(t_) = 0; (Psi)6(t_) = 0;
(Theta)7(t_) = 0; (Phi)7(t_) = 0; (Psi)7(t_) = 0;
(Theta)8(t_) = 0; (Phi)8(t_) = 0; (Psi)8(t_) = 0;
(Theta)9(t_) = 0; (Phi)9(t_) = 0; (Psi)9(t_) = 0;
(Theta)10(t_) = 0; (Phi)10(t_) = 0; (Psi)10(t_) = 0;
(Theta)11(t_) = 0; (Phi)11(t_) = 0; (Psi)11(t_) = 0;
(Theta)12(t_) = 0; (Phi)12(t_) = 0; (Psi)12(t_) = 0;
(Theta)13(t_) = 0; (Phi)13(t_) = 0; (Psi)13(t_) = 0;
(Theta)14(t_) = 0; (Phi)14(t_) = 0; (Psi)14(t_) = 0;
(Theta)15(t_) = 0; (Phi)15(t_) = 0; (Psi)15(t_) = 0;
(Theta)16(t_) = 0; (Phi)16(t_) = 0; (Psi)16(t_) = 0;
(Theta)17(t_) = 0; (Phi)17(t_) = 0; (Psi)17(t_) = 0;
(Theta)18(t_) = 0; (Phi)18(t_) = 0; (Psi)18(t_) = 0;
(Theta)19(t_) = 0; (Phi)19(t_) = 0; (Psi)19(t_) = 0;
(*(Lambda) lengths in cm to the center of mass, measured from
previous node; data from rough measurement on researcher, guessing
center of mass position from DOT paper diagrams*)
(Beta)2 = 0.50044;
(Beta)11 = 0.88494;
(Lambda)2 = 7.2;
(Lambda)3 = 2.6;
(Lambda)4 = 2.6;
(Lambda)5 = 17;
(Lambda)6 = 17;
(Lambda)7 = 11.75;
(Lambda)8 = 11.75;
(Lambda)9 = 7.5;
(Lambda)10 = 7.5;
(Lambda)11 = 6.5;
(Lambda)12 = 19;
(Lambda)13 = 19;
(Lambda)14 = 14;
(Lambda)15 = 14;
(Lambda)16 = 10.5;
(Lambda)17 = 10.5;
(Lambda)18 = 2;
(Lambda)19 = 9;
(*(CapitalLambda)n is the physical length of body segment (n-1) in
cm of the body part that connects to node n, as measured from
previous node; data from rough measurement on researcher *)
(CapitalLambda)2 = 12.5;
(CapitalLambda)3 = 17.25;
(CapitalLambda)4 = 17.25;
(CapitalLambda)5 = 9.7;
(CapitalLambda)6 = 9.7;
(CapitalLambda)7 = 29;
(CapitalLambda)8 = 29;
(CapitalLambda)9 = 27.5;
(CapitalLambda)10 = 27.5;
(CapitalLambda)11 = 11;
(CapitalLambda)12 = 15;
(CapitalLambda)13 = 15;
(CapitalLambda)14 = 41;
(CapitalLambda)15 = 41;
(CapitalLambda)16 = 40;
(CapitalLambda)17 = 40;
(CapitalLambda)18 = 16;
(CapitalLambda)19 = 11;
(* masses of each body segment, from DOT paper *)(*assume abdomen is
fraction of whole torso *)(*assume chest is fraction of whole torso *)
(*assume hips are fraction of whole torso *)
M = {30631*(13/28), 30631*(12/28), 0, 0, 1887, 1794, 1002, 971, 324,
383, 30631*(3/28), 5839, 5601, 2288, 2182, 807, 793, 0, 4025};
(* 3 principal moments of inertia for each segment, units g cm^2;
with the body upright and arms by side and using the abdomen
coordinates the first component is horizontal from left to right, the
second is horizontal front to back, the third is vertical up and down*)
(* data source DOT paper, subject #1, total mass 58.7 kg, height 168
cm*)
(*chest included in abdomen*)
(* hips included in abdomen *)
Ā = {{9315, 14436, 2643}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {146, 132,
23}, {136, 126, 20}, {62, 67, 6}, {52, 54, 6}, {5.3, 4.5,
1.6}, {5.7, 5.7, 1.7}, {0, 0, 0}, {964, 942, 132}, {1086, 1034,
171}, {286, 283, 25}, {310, 290, 35}, {29.6, 5.5, 35.7}, {28.8,
5.6, 30.7}, {0, 0, 0}, {144, 181, 207}};
(*Initial Conditions for the Angular Velocity*)
(Omega)x0 = 0;
(Omega)y0 = 0;
(Omega)z0 = 6;
(* E are relative orientation matrices to adjoining piece that is
closer to center of physical body *)
(* Abdomen serves as reference body segment, its inherent orientation
is given by the identity matrix. For each body part the final
orientation is calculated by rotating that body segment first through
(Psi) (like roll), then (Phi) (like pitch), and finally (Theta)
(like yaw); note that each of the single angle rotation matrices
becomes an identity matrix when that angle is zero *)
a(t_) = {{(Theta)1(t), (Phi)1(t), (Psi)1(t)}, {(Theta)2(
t), (Phi)2(t), (Psi)2(t)}, {(Theta)3(t), (Phi)3(t), (Psi)3(
t)}, {(Theta)4(t), (Phi)4(t), (Psi)4(t)}, {(Theta)5(
t), (Phi)5(t), (Psi)5(t)}, {(Theta)6(t), (Phi)6(t), (Psi)6(
t)}, {(Theta)7(t), (Phi)7(t), (Psi)7(t)}, {(Theta)8(
t), (Phi)8(t), (Psi)9(t)}, {(Theta)9(t), (Phi)9(t), (Psi)9(
t)}, {(Theta)10(t), (Phi)10(t), (Psi)10(t)}, {(Theta)11(
t), (Phi)11(t), (Psi)11(t)}, {(Theta)12(t), (Phi)12(
t), (Psi)12(t)}, {(Theta)13(t), (Phi)13(t), (Psi)13(
t)}, {(Theta)14(t), (Phi)14(t), (Psi)14(t)}, {(Theta)15(
t), (Phi)15(t), (Psi)15(t)}, {(Theta)16(t), (Phi)16(
t), (Psi)16(t)}, {(Theta)17(t), (Phi)17(t), (Psi)17(
t)}, {(Theta)18(t), (Phi)18(t), (Psi)18(t)}, {(Theta)19(
t), (Phi)19(t), (Psi)19(t)}};
E(Theta)(t_) =
Table({{Cos(a(t)((s, 1))), -Sin(a(t)((s, 1))),
0}, {Sin(a(t)((s, 1))), Cos(a(t)((s, 1))), 0}, {0, 0, 1}}, {s,
19});
E(Phi)(t_) =
Table({{1, 0, 0}, {0, Cos(a(t)((s, 2))), -Sin(a(t)((s, 2)))}, {0,
Sin(a(t)((s, 2))), Cos(a(t)((s, 2)))}}, {s, 19});
E(Psi)(t_) =
Table({{Cos(a(t)((s, 3))), 0, -Sin(a(t)((s, 3)))}, {0, 1,
0}, {Sin(a(t)((s, 3))), 0, Cos(a(t)((s, 3)))}}, {s, 19});
R(t_) = Table((E(Theta)(t)((s))).(E(Phi)(t)((s))).(E(Psi)(
t)((s))), {s, 19});
(* orientation matrix relative to abdomen for each body segment *)
B1(t_) = R(t)((1));
B2(t_) = B1(t).R(t)((2));
B3(t_) = B2(t).R(t)((3));
B4(t_) = B2(t).R(t)((4));
B5(t_) = B3(t).R(t)((5));
B6(t_) = B4(t).R(t)((6));
B7(t_) = B5(t).R(t)((7));
B8(t_) = B6(t).R(t)((8));
B9(t_) = B7(t).R(t)((9));
B10(t_) = B8(t).R(t)((10));
B11(t_) = B1(t).R(t)((11));
B12(t_) = B11(t).R(t)((12));
B13(t_) = B11(t).R(t)((13));
B14(t_) = B12(t).R(t)((14));
B15(t_) = B13(t).R(t)((15));
B16(t_) = B14(t).R(t)((16));
B17(t_) = B15(t).R(t)((17));
B18(t_) = B2(t).R(t)((18));
B19(t_) = B18(t).R(t)((19));
B(t_) = {B1(t), B2(t), B3(t), B4(t), B5(t), B6(t), B7(t), B7(t),
B8(t), B9(t), B10(t), B11(t), B12(t), B13(t), B14(t), B15(t),
B16(t), B17(t), B18(t), B19(t)};
(* ê are unit vectors *)
ê0 = {1, 0, 0};
ê1 = {0, 1, 0};
ê2 = {0, 0, 1};
(* ê unit vectors from 1 to 2 connection to shoulders*)
ê3 = {-Sin((Beta)2), 0, Cos((Beta)2)};
ê4 = {Sin((Beta)2), 0, Cos((Beta)2)};
(* ê unit vectors for 1 to 11 connection to hips *)
ê12 = {-Sin((Beta)11), 0, -Cos((Beta)11)};
ê13 = {Sin((Beta)11), 0, -Cos((Beta)11)};
(* u vectors that describe center of mass of segment relative to
previous connector; signs depend on how body part is oriented
relative to connecting parts in the resting body *)
u1 = {0, 0, 0};
u2 = (Lambda)2 ê2;
u3 = -(Lambda)3 ê0;
u4 = (Lambda)4 ê0;
u5 = -(Lambda)5 ê2;
u6 = -(Lambda)6 ê2;
u7 = -(Lambda)7 ê2;
u8 = -(Lambda)8 ê2;
u9 = -(Lambda)9 ê2;
u10 = -(Lambda)10 ê2;
u11 = -(Lambda)11 ê2;
u12 = -(Lambda)12 ê2;
u13 = -(Lambda)13 ê2;
u14 = -(Lambda)14 ê2;
u15 = -(Lambda)15 ê2;
u16 = (Lambda)16 ê1;
u17 = (Lambda)17 ê1;
u18 = (Lambda)18 ê2;
u19 = (Lambda)19 ê2;
(* U vectors are time independent, they let you move from node to
node in the resting body *)
U2 = (CapitalLambda)2 ê2;
U3 = (CapitalLambda)3 ê3;
U4 = (CapitalLambda)4 ê4;
U5 = -(CapitalLambda)5 ê0;
U6 = (CapitalLambda)6 ê0;
U7 = -(CapitalLambda)7 ê2;
U8 = -(CapitalLambda)8 ê2;
U9 = -(CapitalLambda)9 ê2;
U10 = -(CapitalLambda)10 ê2;
U11 = -(CapitalLambda)11 ê2;
U12 = (CapitalLambda)12 ê12;
U13 = (CapitalLambda)13 ê13;
U14 = -(CapitalLambda)14 ê2;
U15 = -(CapitalLambda)15 ê2;
U16 = -(CapitalLambda)16 ê2;
U17 = -(CapitalLambda)17 ê2;
U18 = (CapitalLambda)18 ê2;
U19 = (CapitalLambda)19 ê2;
(CapitalGamma)(t_) = {u1, U2 + B2(t).u2, U2 + B2(t).U3 + u3,
U2 + B2(t).U4 + u4, U2 + B2(t).U3 + B3(t).U5 + B5(t).u5,
U2 + B2(t).U4 + B4(t).U6 + B6(t).u6,
U2 + B2(t).U3 + B3(t).U5 + B5(t).U7 + B7(t).u7,
U2 + B2(t).U4 + B4(t).U6 + B6(t).U8 + B8(t).u8,
U2 + B2(t).U3 + B3(t).U5 + B5(t).U7 + B7(t).U9 + B9(t).u9,
U2 + B2(t).U4 + B4(t).U6 + B6(t).U8 + B8(t).U10 + B10(t).u10,
U11 + B11(t).u11, U11 + B11(t).U12 + B12(t).u12,
U11 + B11(t).U13 + B13(t).u13,
U11 + B11(t).U12 + B12(t).U14 + B14(t).u14,
U11 + B11(t).U13 + B13(t).U15 + B15(t).u15,
U11 + B11(t).U12 + B12(t).U14 + B14(t).U16 + B16(t).u16,
U11 + B11(t).U13 + B13(t).U15 + B15(t).U17 + B17(t).u17,
U11 + B11(t).U13 + B13(t).U15 + B15(t).U17 + B17(t).U18 +
B18(t).u18, U2 + B2(t).U18 + B18(t).U19 + B19(t).u19};
Print("Body Center of Gravity Coordinates=");
(ScriptCapitalC)(
t_) = (Sum((M((s))) (CapitalGamma)(t)((s)), {s,
19}))/(Sum((M((s))), {s, 19}));
(CapitalOmega)(t_) = (CapitalGamma)(t) - (CapitalOmega)(t);
(ScriptCapitalI)(t_) =
Sum(B(t)((s)).DiagonalMatrix({Ā((s, 1)), Ā((s, 2)),
Ā((s, 3))}).Inverse(B(t)((s))) +
M((s)) ((Norm((CapitalOmega)(t)((s)))^2) IdentityMatrix(3) -
Outer(Times, (CapitalOmega)(t)((s)), (CapitalOmega)(
t)((s)))), {s, 19});
Print("Whole-Body Inertia Matrix:((ScriptCapitalI)(t))=");
MatrixForm((ScriptCapitalI)(t));
Print("Eigenvectors((ScriptCapitalI)(t)), in Matrix Form=");
MatrixForm(Eigenvectors((ScriptCapitalI)(t)));
Print("Eigenvalues((ScriptCapitalI)(t))=");
Eigenvalues((ScriptCapitalI)(t));
(*strictly speaking, (Omega)s's are just the list-wise cross product
matricies of (Omega)(t)*)
(Omega)p(s_, t_) = B'(t)((s)).Inverse(B(t)((s)));
(Omega)s(t_) =
Table({(Omega)p(s, t_)((s, 3, 2)), (Omega)p(s, t_)((s, 1,
3)), (Omega)p(s, t_)((s, 2, 1))}, {s, 19});
Lis(t_) =
Table(Cross((CapitalOmega)(t)((s)),
M((s))*(CapitalOmega)'(t)((s))) +
B(t)((s)).DiagonalMatrix({Ā((s, 1)), Ā((s, 2)),
Ā((s, 3))}).Inverse(B(t)((s))).(Omega)s(t)((s)), {s, 19});
Print("Internal/observed Angular Momentum=");
(*The 'N' in 'LiN' stands for 'Net'*)
Li(t_) = Sum(Lis(t)((s)), {s, 19});
Print("Differential Equation Solution =");
sol = NDSolve({(Li'(
t) + (ScriptCapitalI)'(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(
t)} + (ScriptCapitalI)(
t).{(Omega)x'(t), (Omega)y'(t), (Omega)z'(t)} +
Cross({(Omega)x(t), (Omega)y(t), (Omega)z(t)},
Li(t) + (ScriptCapitalI)(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(t)})).ê0 ==
0, (Li'(t) + (ScriptCapitalI)'(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(
t)} + (ScriptCapitalI)(
t).{(Omega)x'(t), (Omega)y'(t), (Omega)z'(t)} +
Cross({(Omega)x(t), (Omega)y(t), (Omega)z(t)},
Li(t) + (ScriptCapitalI)(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(t)})).ê1 ==
0, (Li'(t) + (ScriptCapitalI)'(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(
t)} + (ScriptCapitalI)(
t).{(Omega)x'(t), (Omega)y'(t), (Omega)z'(t)} +
Cross({(Omega)x(t), (Omega)y(t), (Omega)z(t)},
Li(t) + (ScriptCapitalI)(
t).{(Omega)x(t), (Omega)y(t), (Omega)z(t)})).ê2 ==
0, (Omega)x(0) == (Omega)x0, (Omega)y(
0) == (Omega)y0, (Omega)z(0) == (Omega)z0}, {(Omega)x(
t), (Omega)y(t), (Omega)z(t)}, {t, 1, 2},
Method -> {"Automatic", {"EquationSimplification" -> "Residual"}});
VectorQ(s)
(Omega)(t_) = {(Omega)x(t) /.
FilterRules(sol((1)), (Omega)x(t)), (Omega)y(t) /.
FilterRules(sol((1)), (Omega)y(t)), (Omega)z(t) /.
FilterRules(sol((1)), (Omega)z(t)) };
Print("Absolute angular momentum vector, as viewed from the skater's
frame of reference=");
L(t_) = Li(t) + (ScriptCapitalI)(t).(Omega)(t);
(* draw 'frame man' with line segments between nodes, plus extra
lines to shoulders and hips from edges of abdomen region*)
(*Manipulate(Graphics3D({Line({{0,0,0},U2}),Line({U2,U2+B2(t).U3}),
Line({U2+B2(t).U3,U2+B2(t).U3+B3(t).U5}),Line({U2+B2(t).U3+B3(t).U5,
U2+B2(t).U3+B3(t).U5+B5(t).U7}),Line({U2+B2(t).U3+B3(t).U5+B5(t).U7,
U2+B2(t).U3+B3(t).U5+B5(t).U7+B7(t).U9}),Line({U2+B2(t).U3+B3(t).U5+
B5(t).U7+B7(t).U9,U2+B2(t).U3+B3(t).U5+B5(t).U7+B7(t).U9+B9(t).u9}),
Line({U2,U2+B2(t).U4}),Line({U2+B2(t).U4,U2+B2(t).U4+B4(t).U6}),Line({
U2+B2(t).U4+B4(t).U6,U2+B2(t).U4+B4(t).U6+B6(t).U8}),Line({U2+B2(t).
U4+B4(t).U6+B6(t).U8,U2+B2(t).U4+B4(t).U6+B6(t).U8+B8(t).U10}),Line({
U2+B2(t).U4+B4(t).U6+B6(t).U8+B8(t).U10,U2+B2(t).U4+B4(t).U6+B6(t).U8+
B8(t).U10+B10(t).u10}),Line({U2+B2(t).U3,U2+B2(t).U18}),Line({U2+B2(t)
.U4,U2+B2(t).U18}),Line({U2+B2(t).U18,U2+B2(t).U18+B18(t).U19}),Line({
U2+B2(t).U18+B18(t).U19,U2+B2(t).U18+B18(t).U19+B19(t).u19}),Line({{0,
0,0},U11}),Line({U11,U11+B11(t).U12}),Line({U11,U11+B11(t).U13}),Line(
{U11+B11(t).U12,U11+B11(t).U12+B12(t).U14}),Line({U11+B11(t).U13,U11+
B11(t).U13+B13(t).U15}),Line({U11+B11(t).U12+B12(t).U14,U11+B11(t).
U12+B12(t).U14+B14(t).U16}),Line({U11+B11(t).U13+B13(t).U15,U11+B11(t)
.U13+B13(t).U15+B15(t).U17}),Line({U11+B11(t).U12+B12(t).U14+B14(t).
U16,U11+B11(t).U12+B12(t).U14+B14(t).U16+B16(t).u16}),Line({U11+B11(t)
.U13+B13(t).U15+B15(t).U17,U11+B11(t).U13+B13(t).U15+B15(t).U17+B17(t)
.u17}),Line({U11+B11(t).U12,U11+B11(t).u11}),Line({U11+B11(t).U13,U11+
B11(t).u11}),Line({U11,U11+B11(t).u11}),Point((CapitalOmega)(t)),
Blue,Thick,Line({(CapitalOmega)(t),(CapitalOmega)(t)+(10^(-6))(
Eigenvalues((ScriptCapitalI)(t)).ê0)Eigenvectors((ScriptCapitalI)(t)
).ê0}),Green,Thick,Line({(CapitalOmega)(t),(CapitalOmega)(t)+(10^(-
6))(Eigenvalues((ScriptCapitalI)(t)).ê1)Eigenvectors(
(ScriptCapitalI)(t)).ê1}),Red,Thick,Line({(CapitalOmega)(t),
(CapitalOmega)(t)+(10^(-6))(Eigenvalues((ScriptCapitalI)(t)).ê2)
Eigenvectors((ScriptCapitalI)(t)).ê2}),Blue,Thick,Line({
(CapitalOmega)(t),(CapitalOmega)(t)-(10^(-6))(Eigenvalues(
(ScriptCapitalI)(t)).ê0)Eigenvectors((ScriptCapitalI)(t)).ê0}),
Green,Thick,Line({(CapitalOmega)(t),(CapitalOmega)(t)-(10^(-6))(
Eigenvalues((ScriptCapitalI)(t)).ê1)Eigenvectors((ScriptCapitalI)(t)
).ê1}),Red,Thick,Line({(CapitalOmega)(t),(CapitalOmega)(t)-(10^(-6))
(Eigenvalues((ScriptCapitalI)(t)).ê2)Eigenvectors((ScriptCapitalI)(
t)).ê2}),Yellow,Thick,Arrow({(CapitalOmega)(t),(CapitalOmega)(t)+(
10^(-5))Li(t)}),Purple,Thick,Arrow({(CapitalOmega)(t),(CapitalOmega)
(t)+(10^(-6))L(t)}),Orange,Thick,Arrow({(CapitalOmega)(t),
(CapitalOmega)(t)+15(Omega)(t)}),Arrow({(CapitalOmega)13(t),
(CapitalOmega)13(t)+150 (Omega)13(t)})}
,PlotRange(Rule){{-150,150},{-150,150},{-150,150}},SphericalRegion
(Rule)True,RotationAction(Rule)"Clip",AxesEdge(Rule){{-1,-1},{-1,-
1},{-1,-1}},Axes(Rule)True,AxesLabel(Rule){x,y,z}),{t,0,3})*)
```

Oh also, the run time is massive (on my computer anyways) – probably because I’m not yet sure what I’m doing and whilst creating it routinely copied and pasted almost the exact same formula 19 times in many cases, then went in and slightly changed them all by hand. I’ve hit command period each time due to this, not sure if it’s causing my errors but I don’t think so.

So anyways, it says:

a) that some function value is not a list of numbers with dimension 3 at the given point. But I’m stumped as to why, since the output does appear to be three components… unless the <<1>> means that there’s an expression that it hasn’t evaluated, due to an abort, thus, classifying it as not a valid triplet of numbers? I’ve gotten conflicting messages about this point though. The error message just took me back to the Mathematica page for ND Solve, which doesn’t appear to help.

b) I don’t know what Function Value means in this context; I’ve read the documentation, but still only have a fuzzy idea of what it is.

c) Something about *Unable to find initial conditions that satisfy the residual function with specified tolerences.* So I don’t know why it’s saying that, but I *think* it’s secondary to my problem, and I may be able to change that later – but I have a feeling I might not know until I can handle the error message about the .

One other thing, I’ve already read and tried some of the different dif. eq. solving techniques that I managed to find here (In fact if you go to the beginning of the big long code and change (Phi)13(t) from t to 0, you can actually get this as a special case just fine), and just blindly followed some of the error messages so that I have now gotten rid of them. I say this because it reminds me a that, while this is just a conceptual question, why does the code work fine if the guy isn’t moving, but as soon as I make something move it fails to work? I look to forwards to getting better at learning how to use mathematica.

Please tell me if I ought to include more information!

Thanks in advance

**Edit/Update:** So since I started writing this question, I have actually solved the problem of having a massive run time, but now I have completely different problem. It says “Recursion depth of 1024 exceed during the evaluation” (?!? I don’t even think I added anything that was recursive…). (Tell me if I should change the title in light of this)