Can any someone check the following code, I want to plot the streamlines and also I need the figure to be half-spheres only?

lambda=1.5;

eta=1;

beta1=1;beta2=1;

alpha=.01;

eta1=0.1;w=0.1;k=(1/eta.^2).^(0.5); % w=w2/w1

alpha1=((k.^2+(k.^4-4.*alpha.^2).^(0.5))./2).^(0.5);

alpha2=((k.^2-(k.^4-4.*alpha.^2).^(0.5))./2).^(0.5);

a11=(-(alpha1 .^ 2 .* eta - 2 .* eta - 2 .* eta1 - 1) .* alpha1 ./ beta1 .* besselk(0.1e1 ./ 0.2e1, alpha1) - ((alpha1 .^ 2 .* eta) - beta1 - (6 .* eta) - (6 .* eta1) - 0.3e1) .* besselk(0.3e1 ./ 0.2e1, alpha1) ./ beta1) ; ...

a12=(-alpha2 .* (alpha2 .^ 2 .* eta - 2 .* eta - 2 .* eta1 - 1) ./ beta1 .* besselk(0.1e1 ./ 0.2e1, alpha2) - ((alpha2 .^ 2 .* eta) - beta1 - (6 .* eta) - (6 .* eta1) - 0.3e1) .* besselk(0.3e1 ./ 0.2e1, alpha2) ./ beta1) ;

a13= ((alpha1 .^ 2 .* eta - 2 .* eta - 2 .* eta1 - 1) .* alpha1 ./ beta1 .* besseli(0.1e1 ./ 0.2e1, alpha1) - ((alpha1 .^ 2 .* eta) - beta1 - (6 .* eta) - (6 .* eta1) - 0.3e1) .* besseli(0.3e1 ./ 0.2e1, alpha1) ./ beta1) ;

a14=(alpha2 .* (alpha2 .^ 2 .* eta - 2 .* eta - 2 .* eta1 - 1) ./ beta1 .* besseli(0.1e1 ./ 0.2e1, alpha2) - ((alpha2 .^ 2 .* eta) - beta1 - (6 .* eta) - (6 .* eta1) - 0.3e1) .* besseli(0.3e1 ./ 0.2e1, alpha2) ./ beta1) ;

a21=(-0.2e1 .* (eta + eta1) .* alpha1 .* besselk(0.1e1 ./ 0.2e1, alpha1) + (-0.2e1 .* alpha1 .^ 2 .* eta - 0.6e1 .* eta - 0.6e1 .* eta1) .* besselk(0.3e1 ./ 0.2e1, alpha1)) ;

a22=(-0.2e1 .* (eta + eta1) .* alpha2 .* besselk(0.1e1 ./ 0.2e1, alpha2) + (-0.2e1 .* alpha2 .^ 2 .* eta - 0.6e1 .* eta - 0.6e1 .* eta1) .* besselk(0.3e1 ./ 0.2e1, alpha2)) ;

a23=(0.2e1 .* (eta + eta1) .* alpha1 .* besseli(0.1e1 ./ 0.2e1, alpha1) + (-0.2e1 .* alpha1 .^ 2 .* eta - 0.6e1 .* eta - 0.6e1 .* eta1) .* besseli(0.3e1 ./ 0.2e1, alpha1));

a24=(0.2e1 .* (eta + eta1) .* besseli(0.1e1 ./ 0.2e1, alpha2) .* alpha2 + (-0.2e1 .* alpha2 .^ 2 .* eta - 0.6e1 .* eta - 0.6e1 .* eta1) .* besseli(0.3e1 ./ 0.2e1, alpha2));

a31=((alpha1 .^ 2 .* eta .* lambda .^ 2 - lambda .^ 2 - 2 .* eta - 2 .* eta1) .* alpha1 .* (lambda .^ (-0.5e1 ./ 0.2e1)) ./ beta2 .* besselk(0.1e1 ./ 0.2e1, (lambda .* alpha1)) + ((alpha1 .^ 2 .* eta .* lambda .^ 2) + beta2 .* (lambda .^ 3) - (3 .* lambda .^ 2) - (6 .* eta) - (6 .* eta1)) .* (lambda .^ (-0.7e1 ./ 0.2e1)) ./ beta2 .* besselk(0.3e1 ./ 0.2e1, (lambda .* alpha1))) ;

a32= (alpha2 .* (alpha2 .^ 2 .* eta .* lambda .^ 2 - lambda .^ 2 - 2 .* eta - 2 .* eta1) .* (lambda .^ (-0.5e1 ./ 0.2e1)) ./ beta2 .* besselk(0.1e1 ./ 0.2e1, (lambda .* alpha2)) + ((alpha2 .^ 2 .* eta .* lambda .^ 2) + beta2 .* (lambda .^ 3) - (3 .* lambda .^ 2) - (6 .* eta) - (6 .* eta1)) .* (lambda .^ (-0.7e1 ./ 0.2e1)) ./ beta2 .* besselk(0.3e1 ./ 0.2e1, (lambda .* alpha2))) ;

a33= (-(alpha1 .^ 2 .* eta .* lambda .^ 2 - lambda .^ 2 - 2 .* eta - 2 .* eta1) .* alpha1 .* (lambda .^ (-0.5e1 ./ 0.2e1)) ./ beta2 .* besseli(0.1e1 ./ 0.2e1, (lambda .* alpha1)) + ((alpha1 .^ 2 .* eta .* lambda .^ 2) + beta2 .* (lambda .^ 3) - (3 .* lambda .^ 2) - (6 .* eta) - (6 .* eta1)) .* (lambda .^ (-0.7e1 ./ 0.2e1)) ./ beta2 .* besseli(0.3e1 ./ 0.2e1, (lambda .* alpha1))) ;

a34= (-alpha2 .* (alpha2 .^ 2 .* eta .* lambda .^ 2 - lambda .^ 2 - 2 .* eta - 2 .* eta1) .* (lambda .^ (-0.5e1 ./ 0.2e1)) ./ beta2 .* besseli(0.1e1 ./ 0.2e1, (lambda .* alpha2)) + ((alpha2 .^ 2 .* eta .* lambda .^ 2) + beta2 .* (lambda .^ 3) - (3 .* lambda .^ 2) - (6 .* eta) - (6 .* eta1)) .* (lambda .^ (-0.7e1 ./ 0.2e1)) ./ beta2 .* besseli(0.3e1 ./ 0.2e1, (lambda .* alpha2)));

a41=(-0.2e1 .* lambda .^ (-0.3e1 ./ 0.2e1) .* (eta + eta1) .* alpha1 .* besselk(0.1e1 ./ 0.2e1, lambda .* alpha1) - 0.2e1 .* (alpha1 .^ 2 .* eta .* lambda .^ 2 + 0.3e1 .* eta + 0.3e1 .* eta1) .* lambda .^ (-0.5e1 ./ 0.2e1) .* besselk(0.3e1 ./ 0.2e1, lambda .* alpha1)) ;

a42=(-0.2e1 .* lambda .^ (-0.3e1 ./ 0.2e1) .* (eta + eta1) .* alpha2 .* besselk(0.1e1 ./ 0.2e1, lambda .* alpha2) - 0.2e1 .* (alpha2 .^ 2 .* eta .* lambda .^ 2 + 0.3e1 .* eta + 0.3e1 .* eta1) .* lambda .^ (-0.5e1 ./ 0.2e1) .* besselk(0.3e1 ./ 0.2e1, lambda .* alpha2)) ;

a43= (0.2e1 .* lambda .^ (-0.3e1 ./ 0.2e1) .* (eta + eta1) .* alpha1 .* besseli(0.1e1 ./ 0.2e1, lambda .* alpha1) - 0.2e1 .* (alpha1 .^ 2 .* eta .* lambda .^ 2 + 0.3e1 .* eta + 0.3e1 .* eta1) .* lambda .^ (-0.5e1 ./ 0.2e1) .* besseli(0.3e1 ./ 0.2e1, lambda .* alpha1)) ;

a44= (0.2e1 .* lambda .^ (-0.3e1 ./ 0.2e1) .* (eta + eta1) .* alpha2 .* besseli(0.1e1 ./ 0.2e1, lambda .* alpha2) - 0.2e1 .* (alpha2 .^ 2 .* eta .* lambda .^ 2 + 0.3e1 .* eta + 0.3e1 .* eta1) .* lambda .^ (-0.5e1 ./ 0.2e1) .* besseli(0.3e1 ./ 0.2e1, lambda .* alpha2)) ;

A1 = (((a14 .* lambda .* w - a34) .* a43 - a44 .* (a13 .* lambda .* w - a33)) .* a22 + ((-a14 .* lambda .* w + a34) .* a42 + a44 .* (a12 .* lambda .* w - a32)) .* a23 + ((a13 .* lambda .* w - a33) .* a42 - a43 .* (a12 .* lambda .* w - a32)) .* a24) ./ (((-a11 .* a34 + a14 .* a31) .* a43 + (a11 .* a33 - a13 .* a31) .* a44 - a41 .* (-a13 .* a34 + a14 .* a33)) .* a22 + ((a11 .* a34 - a14 .* a31) .* a42 + (-a11 .* a32 + a12 .* a31) .* a44 + a41 .* (-a12 .* a34 + a14 .* a32)) .* a23 + ((-a11 .* a33 + a13 .* a31) .* a42 + (a11 .* a32 - a12 .* a31) .* a43 - a41 .* (-a12 .* a33 + a13 .* a32)) .* a24 - ((a13 .* a34 - a14 .* a33) .* a42 + (-a12 .* a34 + a14 .* a32) .* a43 - a44 .* (-a12 .* a33 + a13 .* a32)) .* a21) ;

B1 = (((-a14 .* lambda .* w + a34) .* a43 + a44 .* (a13 .* lambda .* w - a33)) .* a21 + ((a14 .* lambda .* w - a34) .* a41 - a44 .* (a11 .* lambda .* w - a31)) .* a23 - ((a13 .* lambda .* w - a33) .* a41 - a43 .* (a11 .* lambda .* w - a31)) .* a24) ./ (((a12 .* a34 - a14 .* a32) .* a43 + a44 .* (-a12 .* a33 + a13 .* a32) + a42 .* (-a13 .* a34 + a14 .* a33)) .* a21 + (a41 .* (-a12 .* a34 + a14 .* a32) + (-a11 .* a32 + a12 .* a31) .* a44 - a42 .* (-a11 .* a34 + a14 .* a31)) .* a23 + ((a12 .* a33 - a13 .* a32) .* a41 + (a11 .* a32 - a12 .* a31) .* a43 + (-a11 .* a33 + a13 .* a31) .* a42) .* a24 + a22 .* ((a13 .* a34 - a14 .* a33) .* a41 + (-a11 .* a34 + a14 .* a31) .* a43 - a44 .* (-a11 .* a33 + a13 .* a31))) ;

C1 = (((a14 .* lambda .* w - a34) .* a42 - a44 .* (a12 .* lambda .* w - a32)) .* a21 + ((-a14 .* lambda .* w + a34) .* a41 + a44 .* (a11 .* lambda .* w - a31)) .* a22 + ((a12 .* lambda .* w - a32) .* a41 - a42 .* (a11 .* lambda .* w - a31)) .* a24) ./ ((a42 .* (-a13 .* a34 + a14 .* a33) + a44 .* (-a12 .* a33 + a13 .* a32) - (-a12 .* a34 + a14 .* a32) .* a43) .* a21 + ((a13 .* a34 - a14 .* a33) .* a41 + (a11 .* a33 - a13 .* a31) .* a44 + (-a11 .* a34 + a14 .* a31) .* a43) .* a22 + ((a12 .* a33 - a13 .* a32) .* a41 + (-a11 .* a33 + a13 .* a31) .* a42 - a43 .* (-a11 .* a32 + a12 .* a31)) .* a24 - ((a12 .* a34 - a14 .* a32) .* a41 + a42 .* (-a11 .* a34 + a14 .* a31) - (-a11 .* a32 + a12 .* a31) .* a44) .* a23) ;

D1 = (((-a13 .* lambda .* w + a33) .* a42 + a43 .* (a12 .* lambda .* w - a32)) .* a21 + ((a13 .* lambda .* w - a33) .* a41 - a43 .* (a11 .* lambda .* w - a31)) .* a22 - ((a12 .* lambda .* w - a32) .* a41 - a42 .* (a11 .* lambda .* w - a31)) .* a23) ./ (((a12 .* a34 - a14 .* a32) .* a43 + a44 .* (-a12 .* a33 + a13 .* a32) + a42 .* (-a13 .* a34 + a14 .* a33)) .* a21 + a22 .* ((a13 .* a34 - a14 .* a33) .* a41 + (-a11 .* a34 + a14 .* a31) .* a43 - a44 .* (-a11 .* a33 + a13 .* a31)) + ((a11 .* a34 - a14 .* a31) .* a42 + (-a11 .* a32 + a12 .* a31) .* a44 + a41 .* (-a12 .* a34 + a14 .* a32)) .* a23 + ((a12 .* a33 - a13 .* a32) .* a41 + (-a11 .* a33 + a13 .* a31) .* a42 - a43 .* (-a11 .* a32 + a12 .* a31)) .* a24);

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

a = 1.45/lambda ; %RADIUS

epsilon=.9;L=(1-epsilon).^(1/3);

c =-a/L;

b =a/L;

m =a*35; % NUMBER OF INTERVALS

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

[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);

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;

end

end

end

warning off

v1=(A1 .* besselk(0.3e1 ./ 0.2e1, r .* alpha1) + B1 .* besselk(0.3e1 ./ 0.2e1, r .* alpha2) + C1 .* besseli(0.3e1 ./ 0.2e1, r .* alpha1) + D1 .* besseli(0.3e1 ./ 0.2e1, r .* alpha2)) .* r .^ (-0.1e1 ./ 0.2e1) ;

[DH,h2]=contour(x,y,v1,10,'--b');

hold on

[x,y]=meshgrid([-1:0.001:1]);

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

t=atan2(y,x);

for i1=1:length(x);

for k1=1:length(x);

if sqrt(x(i1,k1).^2+y(i1,k1).^2)>1;

r(i1,k1)=0;

end

end

end

hold on

m1=100;

r1=ones(1,m1+1)*a;r2=ones(1,m1+1)./L;

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

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

set(polar(th,r2,'-k'),'LineWidth',1.3);

%axis square

axis on

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