Do you just need to (A) plot the shape or (B) derive an analytical expression for it? (B) does not belong in this forum, I would think. (A) can be easily solved in the spirit of CFD: set up a grid covering the region of interest, sum up the potential caused by each point singularity (source/sink) at each grid point. Once you have a discretized representation of the potential function, you can differentiate it numerically on this grid to get the velocities, and integrate to get the stream function, and then plot the contours to find closed surfaces. This is sort of like a crude panel method ..