I want to fill the streamlines with another color and fill also the sphere

subplot(2,2,1)

Zeta1=10; %Zeta2=0;

k=.000001;U=1;

alpha1=(sqrt(Zeta1.^2./2+Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

alpha2=(sqrt(Zeta1.^2./2-Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

a = 1 ; %RADIUS

L=.4;

c =-a/L;

b =a/L;

m =a*200; % NUMBER OF INTERVALS

[x,y]=meshgrid([c:(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);

warning on

AA = -(a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 3 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 3 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 3 .* a .* alpha1 .^ 2 + 6 .* a .* alpha1 .* alpha2 + 3 .* a .* alpha2 .^ 2 + 3 .* alpha1 + 3 .* alpha2) .* U ./ alpha2 .^ 2 ./ alpha1 .^ 2;

BB = -0.3e1 .* exp(a .* alpha1) .* sqrt(a .* alpha1) .* (a .* alpha2 + 0.1e1) .* sqrt(0.2e1) .* U ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ (alpha1 - alpha2);

CC = 0.3e1 .* sqrt(0.2e1) .* exp(a .* alpha2) .* (a .* alpha1 + 0.1e1) .* sqrt(a .* alpha2) .* U .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ alpha2 ./ (alpha1 - alpha2);

psi=(BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* sqrt(alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + CC .* sqrt(0.2e1) .* sqrt(pi) .* sqrt(alpha1 .* r) .* exp(-alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* r .* alpha2 .* sqrt(alpha2 .* r) + CC .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha2 .* r) .* r .* alpha1 .* sqrt(alpha1 .* r) + 0.2e1 .* AA .* sqrt(r) .* alpha1 .* sqrt(alpha1 .* r) .* alpha2 .* sqrt(alpha2 .* r)) .* r .^ (-0.3e1 ./ 0.2e1) ./ alpha1 .* (alpha1 .* r) .^ (-0.1e1 ./ 0.2e1) ./ alpha2 .* (alpha2 .* r) .^ (-0.1e1 ./ 0.2e1) .* sin(t) .^ 2 ./ 0.4e1;

[DH,h2]=contour(x,y,psi,5,'k'); %,'ShowText','on'

hold on

m1=100;

r1=ones(1,m1+1)*a;

th=[0:2*pi/m1:2*pi];

set(polar(th,r1,'-k'),'LineWidth',1.1);

title('$\kappa=0.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

axis square

axis on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,2)

Zeta1=10;%Zeta2=0;

k=1;U=1;

alpha1=(sqrt(Zeta1.^2./2+Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

alpha2=(sqrt(Zeta1.^2./2-Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

a = 1 ; %RADIUS

L=.4;

c =-a/L;

b =a/L;

m =a*200; % NUMBER OF INTERVALS

[x,y]=meshgrid([c:(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);

warning on

AA = -(a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 3 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 3 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 3 .* a .* alpha1 .^ 2 + 6 .* a .* alpha1 .* alpha2 + 3 .* a .* alpha2 .^ 2 + 3 .* alpha1 + 3 .* alpha2) .* U ./ alpha2 .^ 2 ./ alpha1 .^ 2;

BB = -0.3e1 .* exp(a .* alpha1) .* sqrt(a .* alpha1) .* (a .* alpha2 + 0.1e1) .* sqrt(0.2e1) .* U ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ (alpha1 - alpha2);

CC = 0.3e1 .* sqrt(0.2e1) .* exp(a .* alpha2) .* (a .* alpha1 + 0.1e1) .* sqrt(a .* alpha2) .* U .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ alpha2 ./ (alpha1 - alpha2);

psi=(BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* sqrt(alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + CC .* sqrt(0.2e1) .* sqrt(pi) .* sqrt(alpha1 .* r) .* exp(-alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* r .* alpha2 .* sqrt(alpha2 .* r) + CC .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha2 .* r) .* r .* alpha1 .* sqrt(alpha1 .* r) + 0.2e1 .* AA .* sqrt(r) .* alpha1 .* sqrt(alpha1 .* r) .* alpha2 .* sqrt(alpha2 .* r)) .* r .^ (-0.3e1 ./ 0.2e1) ./ alpha1 .* (alpha1 .* r) .^ (-0.1e1 ./ 0.2e1) ./ alpha2 .* (alpha2 .* r) .^ (-0.1e1 ./ 0.2e1) .* sin(t) .^ 2 ./ 0.4e1;

[DH,h2]=contour(x,y,psi,5,'k'); %,'ShowText','on'

hold on

m1=100;

r1=ones(1,m1+1)*a;

th=[0:2*pi/m1:2*pi];

set(polar(th,r1,'-k'),'LineWidth',1.1);

title('$\kappa=1.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

axis square

axis on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,3)

Zeta1=10;%Zeta2=0;

k=2;U=1;

alpha1=(sqrt(Zeta1.^2./2+Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

alpha2=(sqrt(Zeta1.^2./2-Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

a = 1 ; %RADIUS

L=.4;

c =-a/L;

b =a/L;

m =a*200; % NUMBER OF INTERVALS

[x,y]=meshgrid([c:(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);

warning on

AA = -(a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 3 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 3 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 3 .* a .* alpha1 .^ 2 + 6 .* a .* alpha1 .* alpha2 + 3 .* a .* alpha2 .^ 2 + 3 .* alpha1 + 3 .* alpha2) .* U ./ alpha2 .^ 2 ./ alpha1 .^ 2;

BB = -0.3e1 .* exp(a .* alpha1) .* sqrt(a .* alpha1) .* (a .* alpha2 + 0.1e1) .* sqrt(0.2e1) .* U ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ (alpha1 - alpha2);

CC = 0.3e1 .* sqrt(0.2e1) .* exp(a .* alpha2) .* (a .* alpha1 + 0.1e1) .* sqrt(a .* alpha2) .* U .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ alpha2 ./ (alpha1 - alpha2);

psi=(BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* sqrt(alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + CC .* sqrt(0.2e1) .* sqrt(pi) .* sqrt(alpha1 .* r) .* exp(-alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* r .* alpha2 .* sqrt(alpha2 .* r) + CC .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha2 .* r) .* r .* alpha1 .* sqrt(alpha1 .* r) + 0.2e1 .* AA .* sqrt(r) .* alpha1 .* sqrt(alpha1 .* r) .* alpha2 .* sqrt(alpha2 .* r)) .* r .^ (-0.3e1 ./ 0.2e1) ./ alpha1 .* (alpha1 .* r) .^ (-0.1e1 ./ 0.2e1) ./ alpha2 .* (alpha2 .* r) .^ (-0.1e1 ./ 0.2e1) .* sin(t) .^ 2 ./ 0.4e1;

[DH,h2]=contour(x,y,psi,5,'k'); %,'ShowText','on'

hold on

m1=100;

r1=ones(1,m1+1)*a;

th=[0:2*pi/m1:2*pi];

set(polar(th,r1,'-k'),'LineWidth',1.1);

title('$\kappa=2.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

axis square

axis on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,4)

Zeta1=10;%Zeta2=0;

k=4;U=1;

alpha1=(sqrt(Zeta1.^2./2+Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

alpha2=(sqrt(Zeta1.^2./2-Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));

a = 1 ; %RADIUS

L=.4;

c =-a/L;

b =a/L;

m =a*200; % NUMBER OF INTERVALS

[x,y]=meshgrid([c:(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);

warning on

AA = -(a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 3 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 3 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 3 .* a .* alpha1 .^ 2 + 6 .* a .* alpha1 .* alpha2 + 3 .* a .* alpha2 .^ 2 + 3 .* alpha1 + 3 .* alpha2) .* U ./ alpha2 .^ 2 ./ alpha1 .^ 2;

BB = -0.3e1 .* exp(a .* alpha1) .* sqrt(a .* alpha1) .* (a .* alpha2 + 0.1e1) .* sqrt(0.2e1) .* U ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ (alpha1 - alpha2);

CC = 0.3e1 .* sqrt(0.2e1) .* exp(a .* alpha2) .* (a .* alpha1 + 0.1e1) .* sqrt(a .* alpha2) .* U .* pi .^ (-0.1e1 ./ 0.2e1) .* a .^ (-0.1e1 ./ 0.2e1) ./ alpha2 ./ (alpha1 - alpha2);

psi=(BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* sqrt(alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + CC .* sqrt(0.2e1) .* sqrt(pi) .* sqrt(alpha1 .* r) .* exp(-alpha2 .* r) .* alpha1 .* alpha2 .* r .^ 2 + BB .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* r .* alpha2 .* sqrt(alpha2 .* r) + CC .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha2 .* r) .* r .* alpha1 .* sqrt(alpha1 .* r) + 0.2e1 .* AA .* sqrt(r) .* alpha1 .* sqrt(alpha1 .* r) .* alpha2 .* sqrt(alpha2 .* r)) .* r .^ (-0.3e1 ./ 0.2e1) ./ alpha1 .* (alpha1 .* r) .^ (-0.1e1 ./ 0.2e1) ./ alpha2 .* (alpha2 .* r) .^ (-0.1e1 ./ 0.2e1) .* sin(t) .^ 2 ./ 0.4e1;

[DH,h2]=contour(x,y,psi,5,'k'); %,'ShowText','on'

hold on

m1=100;

r1=ones(1,m1+1)*a;

th=[0:2*pi/m1:2*pi];

set(polar(th,r1,'-k'),'LineWidth',1.1);

title('$\kappa=4.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

axis square

axis on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%