Just after reading wolfies answer I understood I could calculate the final distribution from the very beginning without all the mid-point steps:

M[x_] := M[x] = Evaluate@FullSimplify@ Integrate[ Boole[a + d + Sqrt[(a - d)^2 + 4 b c] <= x], {a, 0, 1}, {b, 0, 1}, {c, 0, 1}, {d, 0, 1}] gives the CDF and

m[x_] := m[x] = Evaluate@FullSimplify@D[M[x], x] gives the PDF that works perfect with my simulation:

This uses directly the approach of an answer to my previous question.

Just after reading wolfies answer I understood I could calculate the final distribution from the very beginning without all the mid-point steps:

M[x_] := M[x] = Evaluate@FullSimplify@ Integrate[ Boole[a + d + Sqrt[(a - d)^2 + 4 b c] <= x], {a, 0, 1}, {b, 0, 1}, {c, 0, 1}, {d, 0, 1}] gives the CDF and

m[x_] := m[x] = Evaluate@FullSimplify@D[M[x], x] gives the PDF that works perfect with my simulation:

This uses directly the approach of an answer to my previous question.