How can the velocity lines appear around the two spheres?

A =[ -4.7107

0.0012

-0.0056

0.0132

-0.0253

0.0435

-0.0689

0.1031

-0.1473

0.2040

-0.2737

0.3607

-0.4647

0.5927

-0.7425

0.9265

-1.1387

1.4014

-1.7014

2.0810

-2.5114

3.0805

-3.7224

4.6475

-5.6872

7.5039

-9.5388

16.4146

-25.4535

14.3236];

B=[ -3.3794

0.0005

-0.0009

0.0006

-0.0003

0.0001

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000];

C=[ 6.8417

-0.0007

0.0007

-0.0003

0.0001

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000];

AA=[ -4.7100

-0.0012

0.0058

-0.0132

0.0253

-0.0435

0.0689

-0.1032

0.1473

-0.2040

0.2737

-0.3608

0.4648

-0.5928

0.7426

-0.9266

1.1388

-1.4016

1.7016

-2.0813

2.5118

-3.0810

3.7230

-4.6482

5.6881

-7.5050

9.5403

-16.4171

25.4574

-14.3258];

BB=[ -3.3789

-0.0005

0.0009

-0.0006

0.0003

-0.0001

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000];

CC=[ 6.8407

0.0007

-0.0008

0.0003

-0.0001

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000];

a = 1 ; %RADIUS

L=.13;

akm=1;gamma=0.3;arh=0.1;

alphaa=sqrt(((2+akm).*(akm./gamma)).^2+arh.^2);

betaa=(2.*akm.*arh.^2./gamma).^(0.25);

alpha1=sqrt((alphaa.^2+sqrt(alphaa.^4-4.*betaa.^4))./2);

alpha2=sqrt((alphaa.^2-sqrt(alphaa.^4-4.*betaa.^4))./2);

dd=5;

c =-a/L;

b =a/L;

m =a*50; % NUMBER OF INTERVALS

[x,y]=meshgrid([c+dd:(b-c)/m:b],[c:(b-c)/m:b]');

[I J]=find(sqrt(x.^2+y.^2)<(a-.1));

if ~isempty(I);

x(I,J) = 0; y(I,J) = 0;

end

r=sqrt(x.^2+y.^2);

t=atan2(y,x);

r2=sqrt(r.^2+dd.^2-2.*r.*dd.*cos(t));

zet=(r.^2-r2.^2-dd.^2)./(2.*r2.*dd);

warning on

psi1=0;

for i=2:3;

Ai=A(i-1);Bi=B(i-1);Ci=C(i-1);AAi=AA(i-1);BBi=BB(i-1);CCi=C(i-1);

psi1=-psi1-(Ai.*r.^(-i-1)+r.^(-3./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(-3./2).*besselk(i-1./2,r.*alpha2).*Ci).*legendreP(i-1,cos(t))-(AAi.*r2.^(-i-1)+r2.^(-3./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*legendreP(i-1,zet);

end

hold on

[DH1,h1]=contour(x,y,psi1,50,'-k');

%axis square;

title('$(a)$ $\ell=0.1,\;\alpha=1.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

%%%%%%%%%%%%%%%% $\frac{\textstyle a_1+a_2}{\textstyle h}=6.0,\;

hold on

t3 = linspace(0,2*pi,1000);

h2=0;

k2=0;

rr2=1;

x2 = rr2*cos(t3)+h2;

y2 = rr2*sin(t3)+k2;

set(plot(x2,y2,'-k'),'LineWidth',1.1);

fill(x2,y2,'w')

%axis square;

hold on

t2 = linspace(0,2*pi,1000);

h=dd;

k=0;

rr=2;

x1 = rr*cos(t2)+h;

y1 = rr*sin(t2)+k;

set(plot(x1,y1,'-k'),'LineWidth',1.1);

fill(x1,y1,'w')

%axis square;

axis off