mu = 1; % Doublet strength
% GRID:
x = -5:.02:5;
y = -5:.02:5;
for m = 1:length(x)
for n = 1:length(y)
xx(m,n) = x(m); yy(m,n) = y(n);
% Velocity potential function:
phi_Doublet(m,n) = mu * x(m)/(x(m)^2+(y(n)+.01)^2);
% Stream function:
psi_Doublet(m,n) = - mu * y(n)/(x(m)^2+(y(n)+.01)^2);
end
end
% Plots
% Doublet at origin of coordinate system:n
figure(4)
contour(xx,yy,psi_Doublet,[-4:0.25:4],'k'),hold on
contour(xx,yy,phi_Doublet,[-20:0.5:20],'r')
legend('streamlines' ,'potential')
title('Doublet at origin')
axis image,hold off

Doublet at origin