How can I appear the streamline around the spheres?

clc
A=[         -0.2
    0.1
    -0.4
    1.0
    -2.1
    3.9
    -6.9
    11.6
    -18.7
    29.2
    -44.0
    64.9
    -93.5
    132.6
    -184.2
    253.6
    -343.1
    462.3
    -613.4
    815.3
    -1068.2
    1413.3
    -1842.6
    2462.0
    -3232.1
    4517.3
    -6127.8
    10894.9
    -17024.9
    9869.9];
B=[       279.2
    34.9
    -96.3
    177.4
    -248.7
    283.6
    -271.0
    223.3
    -160.5
    102.6
    -58.7
    30.5
    -14.4
    6.3
    -2.5
    0.9
    -0.3
    0.1
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0];
C=[         -4.6
    -0.3
    0.3
    -0.2
    0.1
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0];
AA=[         -1.4
    -0.1
    0.4
    -1.1
    3.0
    -7.5
    17.7
    -40.2
    88.1
    -188.2
    390.9
    -797.9
    1594.7
    -3150.4
    6115.4
    -11791.2
    22393.6
    -42435.9
    79327.3
    -148756.4
    275335.2
    -515215.3
    950981.9
    -1800517.1
    3351792.8
    -6647748.5
    12804165.0
    -32341874.0
    71833152.0
    -59216272.0];
BB=[      15518.5
    -406.6
    1241.4
    -2229.4
    3366.8
    -4306.2
    4732.0
    -4611.8
    4004.3
    -3152.0
    2255.2
    -1486.4
    902.0
    -510.4
    268.4
    -132.9
    61.6
    -27.2
    11.3
    -4.5
    1.7
    -0.6
    0.2
    -0.1
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0];
CC=[        -15.6
    0.3
    -0.4
    0.3
    -0.2
    0.1
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0];
%%%%%%%%%%%%%%%%%%%%
a = 1 ;   %RADIUS
L=.22;
etta=0.01;  %0.02;0.05;0.1d0
u2=1; delta=1.5; ettap=.01; alpha=2.d0; %C=a1/a2=0.1
xi=1./sqrt(etta);xi1=1./sqrt(ettap);
alpha1=sqrt((xi.^2+sqrt(xi.^4-4.*xi.^2.*alpha.^2))./sqrt(2));
alpha2=sqrt((xi.^2-sqrt(xi.^4-4.*xi.^2.*alpha.^2))./sqrt(2));
dd=5;
c =-a/L;
b =a/L;
m =a*40; % 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-.0));
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);
%for i1=1:length(x);
%  for k1=1:length(x);
%    if sqrt(x(i1,k1).^2+y(i1,k1).^2)>1./L;
%        r(i1,k1)=0;r2(i1,k1)=0;
%     end
%  end
%end
warning off
psi1=0;
for i=2:7
    Ai=A(i-1);Bi=B(i-1);Ci=C(i-1);AAi=AA(i-1);BBi=BB(i-1);CCi=CC(i-1);
    psi1=psi1+(Ai.*r.^(-i+1)+r.^(1./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(1./2).*besselk(i-1./2,r.*alpha2).*Ci).*gegenbauerC(i,-1./2, cos(t))+(AAi.*r2.^(-i+1)+r2.^(1./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*gegenbauerC(i,-1./2,zet);
end
hold on
%[DH1,h1]=contour(x,y,psi1,25,'-k','LineWidth',1.1);  %,psi2,'--k',psi2,':k'
%[DH1,h1]=contour(x,y,psi1);
%p1=contour(x,y,psi1,[0.3 0.3],'k','LineWidth',1.1);  %,'ShowText','on'
%p2=contour(x,y,psi1,[0.4 0.4],'r','LineWidth',1.1);  
%p3=contour(x,y,psi1,[0.5 0.5],'g','LineWidth',1.1);  
%p4=contour(x,y,psi1,[0.6 0.6],'b','LineWidth',1.1);  
%p5=contour(x,y,psi1,[0.7 0.7],'c','LineWidth',1.1);  
%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);  
%p7=contour(x,y,psi1,[0.9 0.9],'y','LineWidth',1.1);  
p1=contour(x,y,psi1,[-0.001 0.001],'k','LineWidth',1.1);  %,'ShowText','on'
p2=contour(x,y,psi1,[-0.005 .005],'r','LineWidth',1.1);  
%p3=contour(x,y,psi1,[0.1 0.1],'g','LineWidth',1.1);  
%p4=contour(x,y,psi1,[0.4 0.4],'b','LineWidth',1.1);  
%p5=contour(x,y,psi1,[0.6 0.6],'c','LineWidth',1.1);  
%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);  
%%%%%%%%%%%%%%%  $\frac{\textstyle a_1+a_2}{\textstyle h}=6.0,\;
hold on
t3 = linspace(0,2*pi,1000);
h2=0;
k2=0;
rr2=2;
x2 = rr2*cos(t3)+h2;                   
y2 = rr2*sin(t3)+k2;
set(plot(x2,y2,'-k'),'LineWidth',1.1);
fill(x2,y2,'w')
hold on
t2 = linspace(0,2*pi,1000);
h=dd;
k=0;
rr=1;
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('equal') 
axis on
%xticklabels([])
%yticklabels([])
%legend('0.01','0.05','0.1','0.4','0.6','0.8','Location','northwest')
%title('$\frac{\beta_1}{a_1\mu}=\frac{a_1\beta_2}{\mu}=1.0,\;R_{H}=1.0,\;\frac{a_2}{a_1}=2.0$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')
%title('$(a)\;\; R_{H}=1.0,\;\frac{\kappa}{\mu}=4.0$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')
%%%%%%%%%%%%%%%%%%%%
view([90 90])