Generalizing your mathematics (which are all correct), the 2D stream function automatically satisfies continuity for any 2D case where $\partial\rho / \partial t$ is zero. In flows other than these, you must independently confirm that your results satisfy continuity.
Let's suppose that the boundary is the x-axis. So along the boundary, the stream function is constant. So, $$\frac{\partial^2\psi}{\partial x^2}=0$$ And from the no-slip boundary condition, $$\frac{\partial \psi}{\partial y}=0$$ This can be used to establish the 2nd order finite difference approximation to the value of the vorticity at the boundary:$$\omega=\frac{\partial^2 \psi}{\partial x^2}+\frac{\partial^2 \psi}{\partial y^2}=\frac{2(\psi(I,1)-\psi(I,0))}{(\Delta y)^2}=\omega(I,0)$$
This is used for the boundary condition on $\omega$.
ADDENDUM
First of all, the stream function is known (and constant) at the solid boundaries (because the solid boundaries are stream lines). So it doesn't have to be solved for. It is determined up to an arbitrary constant, and can thus be taken to be zero at one of the boundaries. At the other boundary, the stream function is equal to the volumetric throughput rate per unit width of channel (which is typically known). So you don't need to solve for the stream function at the solid boundaries.
If the tangential derivative $\partial \psi/\partial x=0$ at all locations along the boundary, it's second partial with respect to x must also be equal to zero. This, of course, all follows from the fact that $\psi$ is constant at the boundary.
To integrate the vorticity equation, you need a boundary condition on the vorticity (or at least a 2nd order finite difference approximation to a boundary condition). Just because $\partial \psi/\partial y=0$ does not mean the the second partial of $\psi$ with respect to y is equal to zero at the boundary; this would imply that the vorticity at the boundary is equal to zero, which we know is not correct.
The variable I in the relationships refers to the I'th x grid point. So, back to the boundary condition on vorticity: We have shown so far that, at the boundary, $$\omega=\frac{\partial^2 \psi}{\partial y^2}$$ subject to the constraint that $\partial \psi/\partial y=0$. If we represent these two conditions in 2nd order finite difference form, we obtain:
$$\omega(I,0)=\frac{\psi(I,1)-2\psi(I,0)+\psi(I,-1)}{(\Delta y)^2}$$and$$\frac{\psi(I,1)-\psi(I,-1)}{2\Delta y}=0$$If we combine these two finite difference equations, we obtain a 2nd order finite difference approximation to the value of the vorticity at the boundary:
$$\omega(I,0)=\frac{2(\psi(I,1)-\psi(I,0))}{(\Delta y)^2}$$
I've successfully used this approach to solving these equations many times.
Best Answer
The stream function you give is for a parallel flow plus a source-sink pair located at $x=\pm x_0$. This is not the same thing as a doublet. A doublet is defined as the limit of $x_0\rightarrow0$, $m\rightarrow\infty$ of this source-sink pair such that $m\,x_0=\mbox{const.}=K$
Second, there is no simple closed-form expression of the curve $\Psi=0$ for the case of the Rankine oval that I know of. Numerical computation of that contour will give you the Rankine oval, however, (plus the line $y=0$, of course).
Yes, you are correct, this oval is not an ellipse for any finite choice of the parameters. However, in the above limit of $x_0\rightarrow0$, $m\,x_0=K$ the oval becomes a circle.