deqns = {Subscript[m, 1] x1''[
  t] == ([Lambda]1[t]/Subscript[l, 1]) x1[
   t] - ([Lambda]2[t]/Subscript[l, 2]) (x2[t] - x1[t]),
  Subscript[m, 1] y1''[
  t] == ([Lambda]1[t]/Subscript[l, 1]) y1[
   t] - ([Lambda]2[t]/Subscript[l, 2]) (y2[t] - y1[t]) -
  Subscript[m, 1] g,
  Subscript[m, 2] x2''[
  t] == ([Lambda]2[t]/Subscript[l, 2]) (x2[t] - x1[t]),
  Subscript[m, 2] y2''[
  t] == ([Lambda]2[t]/Subscript[l, 2]) (y2[t] - y1[t]) -
  Subscript[m, 2] g};

aeqns = {x1[t]^2 + y1[t]^2 ==
  Subscript[l, 1]^2, (x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 ==
  Subscript[l, 2]^2};

ics = {x1[0] == 1, y1[0] == 0, y1'[0] == 0, x2[0] == 1, y2[0] == -1,
  y2'[0] == 0};

params = {g -> 9.81, Subscript[m, 1] -> 1, Subscript[m, 2] -> 1,
  Subscript[l, 1] -> 1, Subscript[l, 2] -> 1};

soldp = First[
  NDSolve[{deqns, aeqns, ics} /. params, {x1, y1, x2,
  y2, [Lambda]1, [Lambda]2}, {t, 0, 15000},
  Method -> {"IndexReduction" -> {"Pantelides",
    "ConstraintMethod" -> "Projection"}}]];

SmoothDensityHistogram[
  Map[Function[Evaluate[{x2[#], y2[#]} /. soldp]],
   Range[0, 15000, 0.025]], Mesh -> 5,
  PlotRange -> {{-2, 2}, {-2, 0.1}}]