NDSolve not finding the right solution for a mass hanging of a thread

Question

My NDSolve isn't predicting the right trajectory for one particular case of a problem.

System

Consider a mass vertically hanging by an almost inelastic thread. The thread's force acts acting to the potential

springV[x_] := If[x > l, k (x - l)^2/2, If[x < -l, k (-x - l)^2/2, 0]]

where x is the mass' position. Basically if the mass stretches the string beyond a length l, the string pulls back like a Hook spring$^1$ while if it goes limp, it exerts no force at all.

Let's say the string has an unstretched lengthl. When the mass starts to fall under gravity from initial positions greater than or less than l, NDSolve, gives the right trajectory. However for what is IMO the simplest case $-$ starting from rest at l, NDSolve fails. The solution it produces is one of pure free fall instead of the expected oscillations about equilibrium.

Code

(*aliases*)
echo = Echo;
cls = ClearAll;
first = First;
flatAll[x___] := Flatten[x, Infinity];

(*package - using lagrangian*)
Needs@"VariationalMethods`";

Module[
 {springV, eqns, initConds, x0, v0, x, v, L},
 
 (*the spring potential*)
 springV[x_] := 
  If[x > l, k (x - l)^2/2, If[x < -l, k (-x - l)^2/2, 0]];
  
 (*the lagrangian - T-V; V=0 @x1=0 +ve down*)
 L = m1 x1'[t]^2/2 - springV[x1[t]] + m1 g x1[t];
 
 (*the eqn to solve ma=mg-T*)
 eqns = EulerEquations[L, {x1[t]}, t];
 (*inital conditions*)
 initConds = {x1[0] == x10, x1'[0] == v10};
 (*combine them*)
 eqns = flatAll@{eqns, initConds} // echo;
    (*
    sim params: everything is one except 
    gravity(=9.8)  
    inital velocity(=0)
    k=10^4
*)
 vals = {m1 -> 1, k -> 10^4, l -> 1, x10 -> 1., v10 -> 0, g -> 9.8};
 sol = NDSolve[eqns /. vals // echo, {x1[t]}, {t, 0, 100}];
 xsol = {x1[t]} /. first@sol
 ]

Getting $\ldots$

The yielded solution is a parabola instead of the expected$^2$ oscillations about some $x_0^{[3]}$

Unmagnified, the solution is 1.00098 - 0.00098 Cos[100. t]$^4$ (ignore the plotting artifact around the origin in the magnified expected solution)

Even though in the above I have taken an almost inelastic string by taking a large $k$ relative to $m g$, the same problem exists even with k=1.

Questions

1. Why is NDSolve not giving the right answer and how to correct it?

I have (blindly) tried basic NDSolve variants under "Possible Issues" from docs. to no avail.

---

Test cases (in case your results differ)

In case you have come up with your own code, do test to see it matches the behavior at other x10 values

---

Footer

$^1$ For a hook spring the restorative force is $\propto$ $-$extension$^1$ i.e. its linear.

$^2$using DSolve[{m1 x''[t] == m1 g - k (x[t] - l), x[0] == x10, x'[0] == v10} /. vals, x[t], t]

$^3$$x_0=l+mg/k=1.00098$ for the given sim

$^4$ i.e. it oscillates about the equilibrium position

Answer 1

final modification

potential springV[x] isn't unique at x==1. NDSolve has problems to solve with starting value x1[0]==1.

Try option Method -> {Automatic, "DiscontinuityProcessing" -> False}

sol = NDSolve[eqns /. vals // echo,{x1[t]}, {t, 0, 100} , 
Method -> {Automatic,"DiscontinuityProcessing" -> False}];
xsol = {x1[t]} /. first@sol]   
Plot[xsol, {t, 0, 10}, PlotRange ->{-1.1, 1.1}]