How come this doesn't work as I intended? s[{x0_, y0_}, a_, v_] := NDSolveValue[{ x'[t] == vx[t], WhenEvent[x[t] > 0, vx[t] -> -vx[t]], y'[t] == vy[t] - t, WhenEvent[y[t] -vy[t]], x[0] == x0, vx[0] == v Cos@a, y[0] == y0, vy[0] == v Sin@a}, {x[t], y[t]}, {t, 0, 10}, DiscreteVariables -> {vx, vy}] Simulating falling and bouncing of a 2D mass point. With[{r0 = {-5, 7}}, ParametricPlot[Evaluate@s[r0, -.9, 2], {t, 0, 10}, PlotRange -> {{-10, 10}, {-10, 10}}, Epilog -> {Point@r0}]] I can't understand why it doesn't bounce on the x-axis first, and what does happen there anyway? Event x[t]>0 works fine. Following help I started with lowering the order of equations (and so doubling the number of variables). Here I wrote the equations in a more familiar 2nd order; everything works as expected. I put the function below in a tiny loop which bisected the right initial velocity v0 that returns the ball in its starting position. Could ...