g0 = plot(f(r,0.9,0), (r,-1.5,1.5), legend_label=r'$\theta=0$',
thickness=2, linestyle=':', color='red')
g1 = plot(f(r,0.9,pi/4), (r,-1.5,1.5), legend_label=r'$\theta=\pi/4$',
thickness=2, linestyle='-.', color='grey')
g2 = plot(f(r,0.9,pi/3), (r,-1.5,1.5), legend_label=r'$\theta=\pi/3$',
thickness=2, linestyle='--', color='blue')
g3 = plot(f(r,0.9,pi/2), (r,-1.5,1.5), legend_label=r'$\theta=\pi/2$',
thickness=2, color='violet')
graph = g0+g1+g2+g3
graph.axes_labels([r'$r/m$', r'$\rho^2 (r^2+a^2) + 2 a^2 m r \, \sin^2\theta$'])
graph.set_legend_options(loc='upper right')
graph