Hello There, I've trying to solve 3 nonlinear equation of 3 ball chain connected with springs using Hertz model. Simply 3 balls connect with two springs. I was successful in obtaining the differential equations with no problems. I am using steel balls so you might find parameters for material properties such as Elastic modulus, Poisson ratio, radii of the balls...etc. The main problem is that when using the NDSolve function I get that the system is either singular or stiff. I tried many different methods such StiffnessSwithing, and ExplicitEuler...etc but non of them work especially when I try to solve for larger amount of time such t = 10 seconds. When solving for very small time scale such as micro second that works but definitely not correct. My initial conditions are subject to change with no problem. I would really appreciate your help on this. I am pasting the code in here and I am also going to attach it for you in this post. I did comment on my code so you can understand it with no problems.

(*values of masses,Poisson ratios,modulus of elasticity, Radii of the balls*) 
    values = {m1 -> 1, m2 -> 2, m3 -> 3, E1 -> 210 10^9, E2 -> 210 10^9, 
       E3 -> 210 10^9, \[Nu]1 -> 0.3, \[Nu]2 -> 0.3, \[Nu]3 -> 0.3, 
       R1 -> 6.98*0.01/2, R2 -> 6.98*0.01/2, R3 -> 6.98*0.01/2};
    (*Define the displacement vector for each ball*)
    x[t_] := {x1[t], x2[t], x3[t]};
    (*Assumed initial conditions*)
    IC = {x1[0] == 0.0000001, x2[0] == 0, x3[0] == 0, x1'[0] == 1, 
       x2'[0] == 0, x3'[0] == 0};
    k1 = (1 - \[Nu]1^2)/(\[Pi] E1) /. values;
    k2 = (1 - \[Nu]2^2)/(\[Pi] E2) /. values;
    k3 = (1 - \[Nu]3^2)/(\[Pi] E3) /. values;
    (* Contact Stiffness *)
    K12 = 4/(3 \[Pi] (k1 + k2)) Sqrt[(R1 R2)/(R1 + R2)] /. values;
    K23 = 4/(3 \[Pi] (k2 + k3)) Sqrt[(R2 R3)/(R2 + R3)] /. values;
    (*Elastic deformation at contact points between the balls *)
    \[Delta]12 = (x1[t] - x2[t]) + (R1 + R2) /. values;
    \[Delta]23 = (x2[t] - x3[t]) + (R2 + R3) /. values;
    (* Hertz Equations: Contact forces between the balls *)
    F12 = K12 \[Delta]12^1.5;
    F23 = K23 \[Delta]23^1.5;
    (*Differential equations for each ball*)
    EQ1 = m1 x1''[t] == -F12;
    EQ2 = m2 x2''[t] == F12 - F23;
    EQ3 = m3 x3''[t] == F23;
    (*Put them all together *)
    EOM = {m1 x1''[t] == -F12, m2 x2''[t] == F12 - F23, 
       m3 x3''[t] == F23} /. values
    (*Solve the differential equations *)
    eqn = NDSolve[{EOM, IC}, {x1[t], x2[t], x3[t], x1'[t], x2'[t], 
         x3'[t]}, {t, 0, 10}, Method -> "StiffnessSwitching"] // Flatten //
       Chop
    (* Plot displacement and velocity responses of the system *)
    Plot[Evaluate[x1[t] /. eqn], {t, 0, 10}]
    Plot[Evaluate[x2[t] /. eqn], {t, 0, 10}]
    Plot[Evaluate[x3[t] /. eqn], {t, 0, 10}]

Attachments: