I started with the first of the relevant equations, replacing the p with the operator -iħ∇ and expanding the squared term to yield:
H = (-ħ^2 / 2m)∇^2 + (iqħ/m)A·∇ + (q^2 / 2m)A^2 + qV
But since A = (1/2)B x r
(iqħ/m)A·∇ = (iqħ / 2m)(r x ∇)·B = -(q / 2m)L·B = -(qB_0 / 2m)L_z
and A^2 =...