Nice elegant answer, but when I plotted your solution, the plot was coming out in reverse to Matlab. I needed to reverse the signs of the r and p equations. Do you think the signs should be reversed?
trq = NDSolve[{p'[
t] == -(Simplify[D[H[p, r], r],
r + \[Epsilon] > 0] /. {r -> r[t], p -> p[t]}),
r'[t] == (Simplify[D[H[p, r], p],
r + \[Epsilon] > 0] /. {r -> r[t], p -> p[t]}), r[0] == 1.2,
p[0] == .1}, {r, p}, {t, 18}]
Plot[Evaluate[{r[t], p[t]} /. trq], {t, 0, 18}, PlotRange -> All,
AxesOrigin -> {0, 0}]