Hi, I'm trying to solve a system of 2 equations and 2 variables by using NDsolve function. However, the result is not what I expect. I would like to know if I use the function correctly. This is my first question, I hope this is relevant.

f1[theta1_, theta2_, v1_, 
   v2_] := ((2 + m)*omega^2 * Sin[theta1] + 
     m*omega^2 * Sin[theta1 - 2*theta2] + 
     2*m*(l *v2^2 + v1^2 * Cos[theta1 - theta2]) Sin[
       theta1 - theta2])/ (m * Cos[2*theta1 - 2*theta2] - 2 - m);

f2[theta1_, theta2_, v1_, 
   v2_] := (2*((1 + m) * (v1^2 + omega^2 * Cos[theta1]) + 
       l*m*v2^2 * Cos[theta1 - theta2])* 
     Sin[theta1 - theta2]) / (-l*(m * Cos[2*theta1 - 2*theta2]) - 2 - 
     m);
m1 = 0.5;
m2 = 2;
l1 = 0.6;
l2 = 0.9;
omega=0.5;

   solNl = NDSolve[{
   th3''[s] == f1[th3[s], th4[s], th3'[s], th4'[s]],
   th4''[s] == f2[th3[s], th4[s], th3'[s], th4'[s]],
   th3[0] == 0.6,
   th4[0] == 0.6,
   th3'[0] == 0,
   th4'[0] == 0}, {th3, th4}, {s, 0, 30 Pi}];

   Sol3[t_] := th3[t] /. solNl;
   Sol3[8 Pi]

As result I got

NDSolve::ndsz: At s == 24.140269374312506`, step size is effectively zero; singularity or stiff     system suspected.
InterpolatingFunction::dmval: Input value {8 \[Pi]} lies outside the range of data in the interpolating function. Extrapolation will be used.
{3.72634*10^49}

Any help will be appreciated.