How to plot two concentric spheres with radius1=1 and radius2 =10 in MATLAB and filled with fluid flow (I have one velocity component only)

How to plot two concentric spheres with radius1=1 and radius2 =10 in MATLAB and filled with fluid flow (I have one velocity component only)

[Dyuman Joshi]

What have you attempted yet?

What exactly does fluid flow mean?

[shreen elsapa]

I have only the particles velocity, can I plot the streamlines by this function?

velocity = (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) .* sin(theta);

[Dyuman Joshi]

What are the values of the variables?

You will need atleast 2 velocity components to plot streamlines

[shreen elsapa]

%I am try to plot this but the streamlines did not appears between two concentric???

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

a = 1 ; %RADIUS

%for different alpha

lambda=10;

eta=1;

beta1=10.^10;beta2=10.^10;

alpha=.00001;

eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5); a=1; % 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);

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

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

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-.2));

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,12,'-k');

hold on

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

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 off

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

Answer 1

[Torsten]

lambda=1.5;
eta=1;
beta1=10;beta2=10;
alpha=.00001;
eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5);  a=1;   % 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Rmin = 1;
Rmax = 2.5;
range = -Rmax:0.01:Rmax;
[x,y] = meshgrid(range);
r = sqrt(x.^2+y.^2);
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) ;
v1(r<Rmin)=NaN;
v1(r>Rmax)=NaN;
contourf(x,y,v1)
hold on
theta = linspace(0,2*pi,100);
xfill = Rmin*cos(theta);
yfill = Rmin*sin(theta);
fill(xfill,yfill,[0.7 0.7 0.7])
hold on
theta = linspace(0,pi/2,25);
xfill = [0,Rmax,Rmax,Rmax*cos(theta)];
yfill = [Rmax,Rmax,0,Rmax*sin(theta)];
fill(xfill,yfill,[0.7 0.7 0.7])
hold on
theta = linspace(pi/2,pi,25);
xfill = [-Rmax,-Rmax,0,Rmax*cos(theta)];
yfill = [0,Rmax,Rmax,Rmax*sin(theta)];
fill(xfill,yfill,[0.7 0.7 0.7])
hold on
theta = linspace(pi,3/2*pi,25);
xfill = [0,-Rmax,-Rmax,Rmax*cos(theta)];
yfill = [-Rmax,-Rmax,0,Rmax*sin(theta)];
fill(xfill,yfill,[0.7 0.7 0.7])
hold on
theta = linspace(3/2*pi,2*pi,25);
xfill = [Rmax,Rmax,0,Rmax*cos(theta)];
yfill = [0,-Rmax,-Rmax,Rmax*sin(theta)];
fill(xfill,yfill,[0.7 0.7 0.7])
hold off
colorbar

[shreen elsapa]

thanks so much for trying to help me, please see the below I tryed to plot but I need to plot half of both spheres only, what can I do?

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

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

[Torsten]

lambda=1.5;
eta=1;
beta1=10;beta2=10;
alpha=.00001;
eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5);  a=1;   % 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rmin = 1;
Rmax = 2.5;
range = -Rmax:0.01:Rmax;
[x,y] = meshgrid(range);
r = sqrt(x.^2+y.^2);
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) ;
v1(r<Rmin)=NaN;
v1(r>Rmax)=NaN;
v1(y<0) = NaN;
contourf(x,y,v1)
hold on
theta = linspace(0,pi,100);
xfill = Rmin*cos(theta);
yfill = Rmin*sin(theta);
fill([-Rmin,Rmin,xfill],[0,0,yfill],[0.7 0.7 0.7])
hold on
theta = linspace(0,pi/2,25);
xfill = [0,Rmax,Rmax,Rmax*cos(theta)];
yfill = [Rmax,Rmax,0,Rmax*sin(theta)];
fill(xfill,yfill,[0.7 0.7 0.7])
hold on
theta = linspace(pi/2,pi,25);
xfill = [-Rmax,-Rmax,0,Rmax*cos(theta)];
yfill = [0,Rmax,Rmax,Rmax*sin(theta)];
fill(xfill,yfill,[0.7 0.7 0.7])
hold off
colorbar
axis equal
ylim([0 Rmax])

Answer 2

[Torsten]

I'm surprised that velocity increases with growing r. Continuity equation should force velocity to decrease. But maybe my parameters are wrong or I didn't understand what is plotted.

r = 1:0.1:10;
A1 = 1;
B1 = 1;
C1 = 1;
D1 = 1;
alpha1 = 1;
alpha2 = 1;
theta = pi/2;
velocity = @(r)(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) .* sin(theta);
plot(r,velocity(r))

[shreen elsapa]

%I am try to plot this but the streamlines did not appears between two concentric???

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

a = 1 ;   %RADIUS
%for different alpha
lambda=10;
eta=1;
beta1=10.^10;beta2=10.^10;
alpha=.00001;
eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5);  a=1;   % 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon=.9;L=(1-epsilon).^(1/3);
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-.2));
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,12,'-k');
hold on
[x,y]=meshgrid([-1:0.002: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)>a;
            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 off

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

[shreen elsapa]

please see the previous code

[Torsten]

please see the previous code

What do you mean ? That's what I get from your code.

Nothing is plotted since v1 = NaN identically.

[shreen elsapa]

lambda=1.5;

eta=1;

beta1=10;beta2=10;

alpha=.00001;

eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5); a=1; % 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/lambda ; %RADIUS

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

c =-a/L;

b =a/L;

m =a*20; % 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-0.08));

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,50,'--k');

hold on

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

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 off

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

[shreen elsapa]

Thanks. I want the color in the annulus between two spheres, if you can?

[Torsten]

Done.

[shreen elsapa]

thanks, please if I want to remove the gry color from the background

lambda=2.5;
eta=1;
%beta1=10.^10;beta2=10.^10;
beta1=10.^10;beta2=10.^10;
alpha=1;
eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5);  a=1;   % 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rmin = 1;
Rmax = 2.5;
range = -Rmax:0.01:Rmax;
[x,y] = meshgrid(range);
r = sqrt(x.^2+y.^2);
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) ;
v1(r<Rmin)=NaN;
v1(r>Rmax)=NaN;
v1(y<0) = NaN;
contour(x,y,v1,'ShowText','on')
hold on
theta = linspace(0,pi,100);
xfill = Rmin*cos(theta);
yfill = Rmin*sin(theta);
fill([-Rmin,Rmin,xfill],[0,0,yfill],[0.7 0.7 0.7])
hold off
colorbar
axis equal
ylim([0 Rmax])

[Torsten]

Done.

But your velocity profile does not satisfy the continuity equation: v(r)*(4*pi*r^2) = constant

[shreen elsapa]

velocity equal constant means the stream function, is this you mean?

[Torsten]

I mean that the volume flowing through the surfaces of a sphere must remain constant for all radii - assuming that the density remains constant and that there are no mass sources inside the sphere.

With only one velocity component, the volume flow is Vdot = v1(r)*4*pi*r^2, and this value should remain constant.

[shreen elsapa]

Ok. Can I plot the outer contour of the outside sphere solid with black color to be different from the other countours

[Torsten]

Rmax = 2.5;
theta = linspace(0,pi,100);
xsphere = Rmax*cos(theta);
ysphere = Rmax*sin(theta);
plot(xsphere,ysphere,'Color',"k")
axis equal
xlim([-Rmax Rmax])
ylim([0 Rmax])

I think you should be able to include this in the remaining code.

[shreen elsapa]

Dear, you are right the velocity should be decayed with th eincrease of r , there are many physical parameters has the effect on the velocity that did the required effect.

Thanks for your comments,

[Torsten]

Decay is visible in your contours plot, but not the correct one. The following plot should be a horizontal line. But maybe I don't understand what your v1 really is - I assumed it is some kind of velocity in a sphere in radial direction.

lambda=2.5;
eta=1;
%beta1=10.^10;beta2=10.^10;
beta1=10.^10;beta2=10.^10;
alpha=1;
eta1=0.0;w=0.1;k=(1/eta.^2).^(0.5);  a=1;   % 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rmin = 1;
Rmax = 2.5;
r = Rmin:0.01:Rmax;
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) ;
plot(r,4*pi*r.^2.*v1)

[shreen elsapa]

Yes, there are some parameters that did variance because the motion assumed rotational and there are other impacts on this stream function.

Again thanks for your help.