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[{θ''[t] + φ'[t]*Sin[φ[t]]*(Cos[θ[t]]*Cos[t] + 
          Sin[θ[t]]*Sin[t]) + Cos[φ[t]]*(θ'[t]*Sin[θ[t]]*Cos[t] + 
          Cos[θ[t]]*Sin[t] - θ'[t]*Cos[θ[t]]*Sin[t] - Sin[θ[t]]*Cos[t]) == 
          (φ'[t]*Sin[φ[t]]*(Cos[θ[t]]*Cos[t] + Sin[θ[t]]*Sin[t]) - θ'[t]*
          Cos[φ[t]]*(Cos[θ[t]]*Sin[t] - Sin[θ[t]]*Cos[t])) + (φ'[t]^2)*
          Cos[θ[t]]*Sin[θ[t]] - 9.81*Sin[θ[t]], φ'[t]*Cos[φ[t]]*(Sin[θ[t]]*Cos[t] - 
          Cos[θ[t]]*Sin[t]) + Sin[φ[t]]*(θ'[t]*Cos[θ[t]]*Cos[t] - 
          Sin[θ[t]]*Sin[t] + θ'[t]*Sin[θ[t]]*Sin[t] - Cos[θ[t]]*Cos[t]) + 
          φ''[t]*((Sin[θ[t]])^2) + 2*θ'[t]*φ'[t]*Sin[θ[t]]*Cos[θ[t]] == 
          ((φ'[t]^2)*Cos[φ[t]]*(Sin[θ[t]]*Cos[t] - Cos[θ[t]]*Sin[t]) + θ'[t]*φ'[t]*
          Sin[φ[t]]*(Cos[θ[t]]*Cos[t] + Sin[θ[t]]*Sin[t])), 
          θ[0] == Pi/8, φ[0] == Pi/8, θ'[0] == 0, φ'[0] == 0}, {θ, φ}, {t, 0.5}]
    
    Plot[Evaluate[{θ[t], φ[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.

----

Whatever I do, it is possible that my equations could be wrong - they're derived from this Lagrangian:

    l = φ'[t] Sin[φ[t]] (Sin[θ[t]] Cos[t] - Cos[θ[t]] Sin[t]) - θ'[ t] Cos[φ[t]]             
        (Cos[θ[t]] Cos[t] + Sin[θ[t]] Sin[t]) + 0.5*(θ'[t]^2 + φ'[t]^2 Sin[θ[t]]^2 + 1) + 
         9.81*Cos[θ[t]]