Given a divergenceless vector field $\mathbf{B}(r,\theta)$ defined on an uneven 2D numerical grid in $(r,\theta)$ spherical coordinates, I'd like to plot its streamlines.
One way to do it is to plot level curves (contours) of the Stokes stream function $\psi(r,\theta)$, determined by the following equations: $$B_r=\frac{1}{r^2\sin\theta}\frac{\partial\psi}{\partial\theta}\;\;\;\;B_\theta=-\frac{1}{r\sin\theta}\frac{\partial\psi}{\partial r}$$
Since I don't have any experience with numerically solving PDEs, I'm not sure how to approach the problem. I've looked at some fluid mechanics texts with numerical methods but I couldn't find what I need.
Thanks.