I was wondering if there was a way to help the solutions to these differential equations not hit a singularity as quickly as they currently do. The equations are very complex (derived from a Lagrangian), so it could be that Mathematica is struggling because of that. Depending on the initial conditions, I'm only able to solve the equations for values of time (t) up to half a second or a few seconds if lucky. This is my input: s = NDSolve[{\[Theta]''[t] + \[CurlyPhi]'[t]* Sin[\[CurlyPhi][t]]*(Cos[\[Theta][t]]*Cos[t] + Sin[\[Theta][t]]*Sin[t]) + Cos[\[CurlyPhi][t]]*(\[Theta]'[t]*Sin[\[Theta][t]]*Cos[t] + Cos[\[Theta][t]]*Sin[t] - \[Theta]'[t]*Cos[\[Theta][t]]* Sin[t] - Sin[\[Theta][t]]*Cos[t]) == (\[CurlyPhi]'[t]* Sin[\[CurlyPhi][t]]*(Cos[\[Theta][t]]*Cos[t] + Sin[\[Theta][t]]*Sin[t]) - \[Theta]'[t]* Cos[\[CurlyPhi][t]]*(Cos[\[Theta][t]]*Sin[t] - Sin[\[Theta][t]]*Cos[t])) + (\[CurlyPhi]'[t]^2)* Cos[\[Theta][t]]*Sin[\[Theta][t]] - 9.81*Sin[\[Theta][t]], \[CurlyPhi]'[t]* Cos[\[CurlyPhi][t]]*(Sin[\[Theta][t]]*Cos[t] - Cos[\[Theta][t]]*Sin[t]) + Sin[\[CurlyPhi][t]]*(\[Theta]'[t]*Cos[\[Theta][t]]*Cos[t] - Sin[\[Theta][t]]*Sin[t] + \[Theta]'[t]*Sin[\[Theta][t]]* Sin[t] - Cos[\[Theta][t]]*Cos[t]) + \[CurlyPhi]''[ t]*((Sin[\[Theta][t]])^2) + 2*\[Theta]'[t]*\[CurlyPhi]'[t]*Sin[\[Theta][t]]* Cos[\[Theta][t]] == ((\[CurlyPhi]'[t]^2)* Cos[\[CurlyPhi][t]]*(Sin[\[Theta][t]]*Cos[t] - Cos[\[Theta][t]]*Sin[t]) + \[Theta]'[t]*\[CurlyPhi]'[t]* Sin[\[CurlyPhi][t]]*(Cos[\[Theta][t]]*Cos[t] + Sin[\[Theta][t]]*Sin[t])), \[Theta][0] == Pi/8, \[CurlyPhi][0] == Pi/8, \[Theta]'[0] == 0, \[CurlyPhi]'[0] == 0}, {\[Theta], \[CurlyPhi]}, {t, 0.5}] Plot[Evaluate[{\[Theta][t], \[CurlyPhi][t]} /. s], {t, 0, 0.5}, PlotStyle -> Automatic] Sorry if this is a little messy, but I'm simply wanting to know whether Mathematica can do better than come up with a better numerical solution than one which is only reasonable for one second or so.