Correction: evaluation of Matérn 3/2 GP density in Stan using the SDE representation takes quadratic time!
Summary: The method proposed in http://www.juhokokkala.fi/blog/posts/linear-time-evaluation-of-matern-32-gp-density-in-stan-using-the-sde-representation/ has quadratic time complexity (per evaluation of density+gradient) due to the gradient evaluation. This explains why the speedup seen in the experiment was only proportional to the number of input points (consistent with the difference between cubic and quadratic time complexity), not to the square of the number of input points as I expected. NB: I have not examined the code Stan produces, the possibility remains that the time-complexity is even higher due to something I have missed...