Why the streamlines do not appear around the sphere?

clc
a = 1; % radius (m)
b = a*2;
c = -a*2;
n = a*15; % number of intervals
% the spherical coordinate solution for a falling sphere looking only at x-y plane (looking from side not top)
[x,y]=meshgrid((c:(b-c)/n:b),(c:(b-c)/n:b)');
% Preliminar DATA & purification
for i = 1:length(x)
    for k = 1:length(x)
        if sqrt(x(i,k).^2 + y(i,k).^2) < a
            x(i,k) = 0;
            y(i,k) = 0;
        end
    end
end
% definition of polar variables
r = sqrt(x.^2 + y.^2);
theta = atan2(y,x);
% creation of the streamline function
beta1=1; beta2=1;
Zeta1=0.01;Zeta2=0.01;
k=1;U=1;
alpha1=real(sqrt(Zeta1.^2./2+Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));
alpha2=real(sqrt(Zeta1.^2./2-Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));
AA = -(Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 4 .* alpha2 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 4 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .* beta1 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .* alpha2 .^ 2 .* beta1 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 2 .* beta1 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 2 .* alpha2 .^ 3 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 2 .* beta2 + 2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 4 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 4 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .* beta1 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 .* beta2 + 4 .* Zeta1 .^ 4 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 2 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 3 + 20 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .* alpha2 .^ 3 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 4 + 8 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta2 + 12 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 .* beta2 + 8 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 12 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* beta1 + 40 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 36 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 + 36 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 3 .* beta1 + 4 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 4 .* alpha2 .^ 2 + 8 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 4 + 8 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 4 .* alpha2 .^ 2 + 16 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .* alpha1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .* alpha2 .* beta2 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 12 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .^ 2 .* beta1 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .^ 2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 + 72 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta1 + 72 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta1 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 .^ 2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .^ 3 + 24 .* Zeta1 .^ 4 .* a .* alpha1 .^ 2 + 48 .* Zeta1 .^ 4 .* a .* alpha1 .* alpha2 + 12 .* Zeta1 .^ 4 .* a .* alpha1 .* beta1 + 24 .* Zeta1 .^ 4 .* a .* alpha2 .^ 2 + 12 .* Zeta1 .^ 4 .* a .* alpha2 .* beta1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 + 96 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* beta1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .* beta1 + 24 .* Zeta1 .^ 4 .* alpha1 + 24 .* Zeta1 .^ 4 .* alpha2 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* alpha1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* alpha2) .* a ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta2 + 2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 4 .* Zeta1 .^ 4 .* a .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 + 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 8 .* Zeta1 .^ 4 .* a + 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 + 8 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 + 4 .* Zeta1 .^ 2 .* a .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 + 16 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 + 8 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 + 4 .* Zeta1 .^ 2 .* alpha1 + 4 .* Zeta1 .^ 2 .* alpha2 + 8 .* Zeta2 .^ 2 .* alpha1 + 8 .* Zeta2 .^ 2 .* alpha2) ./ alpha2 .^ 2 ./ alpha1 .^ 2;
BB = -0.3e1 .* exp(alpha1 .* a) .* sqrt(alpha1 .* a) .* (Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* beta2 + 0.4e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha2 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha2 + 0.4e1 .* Zeta1 .^ 2 + 0.8e1 .* Zeta2 .^ 2) .* sqrt(a) .* Zeta1 .^ 2 .* (a .* beta1 + 0.2e1) .* sqrt(0.2e1) ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 .* beta2 - Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 .* beta2 + 0.2e1 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 - 0.2e1 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 + 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 - 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 .* beta2 + 0.4e1 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 - 0.4e1 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 + 0.4e1 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* beta1 - 0.4e1 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .* beta1 + 0.12e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 0.8e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 - 0.12e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 - 0.8e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 .* beta2 + 0.4e1 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 0.4e1 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 0.8e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 0.8e1 .* Zeta1 .^ 4 .* a .* alpha1 - 0.8e1 .* Zeta1 .^ 4 .* a .* alpha2 + 0.16e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 - 0.16e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .^ 3 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* a .* alpha2 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 0.8e1 .* Zeta2 .^ 2 .* a .* alpha2 .^ 3 + 0.4e1 .* Zeta1 .^ 2 .* alpha1 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* alpha2 .^ 2 + 0.8e1 .* Zeta2 .^ 2 .* alpha1 .^ 2 - 0.8e1 .* Zeta2 .^ 2 .* alpha2 .^ 2);
CC = 0.3e1 .* (Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 4 .* Zeta1 .^ 2 .* a .* alpha1 + 8 .* Zeta2 .^ 2 .* a .* alpha1 + 4 .* Zeta1 .^ 2 + 8 .* Zeta2 .^ 2) .* sqrt(a) .* sqrt((alpha2 .* a)) .* (Zeta1 .^ 2) .* (a .* beta1 + 2) .* exp((alpha2 .* a)) .* sqrt(0.2e1) ./ alpha2 .* pi .^ (-0.1e1 ./ 0.2e1) ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 .* beta2 - Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 - 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 - 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 - 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 - 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* beta1 - 4 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 - 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 - 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 8 .* Zeta1 .^ 4 .* a .* alpha1 - 8 .* Zeta1 .^ 4 .* a .* alpha2 + 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 - 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 3 + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 4 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 4 .* Zeta1 .^ 2 .* a .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 3 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 8 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 8 .* Zeta2 .^ 2 .* a .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* alpha1 .^ 2 - 4 .* Zeta1 .^ 2 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* alpha1 .^ 2 - 8 .* Zeta2 .^ 2 .* alpha2 .^ 2);
psi = (BB .* sqrt(alpha2 .* r) .* sqrt(0.2e1) .* sqrt(pi) .* exp(-alpha1 .* r) .* alpha1 .* alpha2 .* r .^ 2 + CC .* exp(-alpha2 .* r) .* sqrt(0.2e1) .* sqrt(pi) .* sqrt(alpha1 .* 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(theta) .^ 2 ./ 0.4e1;
% laminar stream function in polar coordinate
% stream line Plot
contour(x,y,psi,25,'Linewidth',1.2);  hold on; axis square;

Warning: Contour not rendered for non-finite ZData

% circle plot
r = ones(1,n+1)*a;
t = (0:2*pi/n:2*pi);
polar(t,r,'k');
r = sqrt(x.^2 + y.^2);
theta = atan2(y,x);
% creating The Filled Circle
t_r  = (0:.01:2*pi);
xxx = a*cos(t_r);
yyy = a*sin(t_r);
fill(xxx,yyy,'y'); axis square

[Mehmed Saad]

the variable psi in line 29 is NaN most probably.

[Shreen El-Sapa]

beta1=0.00001; beta2=0.00001;

Zeta1=1;Zeta2=1;

k=.00001;U=1;

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

alpha2=real(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 = -(Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 4 .* alpha2 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 4 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .* beta1 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .* alpha2 .^ 2 .* beta1 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 2 .* beta1 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 2 .* alpha2 .^ 3 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 2 .* beta2 + 2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 4 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 4 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .* beta1 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 .* beta2 + 4 .* Zeta1 .^ 4 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 2 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 3 + 20 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .* alpha2 .^ 3 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 4 + 8 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta2 + 12 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 .* beta2 + 8 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 12 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* beta1 + 40 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 36 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 + 36 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 3 .* beta1 + 4 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 4 .* alpha2 .^ 2 + 8 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 4 + 8 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 4 .* alpha2 .^ 2 + 16 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .* alpha1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .* alpha2 .* beta2 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 12 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .^ 2 .* beta1 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .^ 2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 + 72 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta1 + 72 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta1 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 .^ 2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .^ 3 + 24 .* Zeta1 .^ 4 .* a .* alpha1 .^ 2 + 48 .* Zeta1 .^ 4 .* a .* alpha1 .* alpha2 + 12 .* Zeta1 .^ 4 .* a .* alpha1 .* beta1 + 24 .* Zeta1 .^ 4 .* a .* alpha2 .^ 2 + 12 .* Zeta1 .^ 4 .* a .* alpha2 .* beta1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 + 96 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* beta1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .* beta1 + 24 .* Zeta1 .^ 4 .* alpha1 + 24 .* Zeta1 .^ 4 .* alpha2 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* alpha1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* alpha2) .* a ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta2 + 2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 4 .* Zeta1 .^ 4 .* a .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 + 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 8 .* Zeta1 .^ 4 .* a + 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 + 8 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 + 4 .* Zeta1 .^ 2 .* a .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 + 16 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 + 8 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 + 4 .* Zeta1 .^ 2 .* alpha1 + 4 .* Zeta1 .^ 2 .* alpha2 + 8 .* Zeta2 .^ 2 .* alpha1 + 8 .* Zeta2 .^ 2 .* alpha2) ./ alpha2 .^ 2 ./ alpha1 .^ 2;

BB = -0.3e1 .* exp(alpha1 .* a) .* sqrt(alpha1 .* a) .* (Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* beta2 + 0.4e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha2 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha2 + 0.4e1 .* Zeta1 .^ 2 + 0.8e1 .* Zeta2 .^ 2) .* sqrt(a) .* Zeta1 .^ 2 .* (a .* beta1 + 0.2e1) .* sqrt(0.2e1) ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 .* beta2 - Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 .* beta2 + 0.2e1 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 - 0.2e1 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 + 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 - 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 .* beta2 + 0.4e1 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 - 0.4e1 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 + 0.4e1 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* beta1 - 0.4e1 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .* beta1 + 0.12e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 0.8e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 - 0.12e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 - 0.8e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 .* beta2 + 0.4e1 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 0.4e1 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 0.8e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 0.8e1 .* Zeta1 .^ 4 .* a .* alpha1 - 0.8e1 .* Zeta1 .^ 4 .* a .* alpha2 + 0.16e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 - 0.16e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .^ 3 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* a .* alpha2 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 0.8e1 .* Zeta2 .^ 2 .* a .* alpha2 .^ 3 + 0.4e1 .* Zeta1 .^ 2 .* alpha1 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* alpha2 .^ 2 + 0.8e1 .* Zeta2 .^ 2 .* alpha1 .^ 2 - 0.8e1 .* Zeta2 .^ 2 .* alpha2 .^ 2);

CC =0.3e1 .* (Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 4 .* Zeta1 .^ 2 .* a .* alpha1 + 8 .* Zeta2 .^ 2 .* a .* alpha1 + 4 .* Zeta1 .^ 2 + 8 .* Zeta2 .^ 2) .* sqrt(a) .* sqrt((alpha2 .* a)) .* (Zeta1 .^ 2) .* (a .* beta1 + 2) .* exp((alpha2 .* a)) .* sqrt(0.2e1) ./ alpha2 .* pi .^ (-0.1e1 ./ 0.2e1) ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 .* beta2 - Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 - 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 - 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 - 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 - 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* beta1 - 4 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 - 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 - 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 8 .* Zeta1 .^ 4 .* a .* alpha1 - 8 .* Zeta1 .^ 4 .* a .* alpha2 + 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 - 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 3 + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 4 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 4 .* Zeta1 .^ 2 .* a .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 3 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 8 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 8 .* Zeta2 .^ 2 .* a .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* alpha1 .^ 2 - 4 .* Zeta1 .^ 2 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* alpha1 .^ 2 - 8 .* Zeta2 .^ 2 .* alpha2 .^ 2);

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

p1=contour(x,y,psi,[-.1 -.1],'k','LineWidth',1.1); %,'ShowText','on'

p2=contour(x,y,psi,[-.3 -.3],'r','LineWidth',1.1);

p3=contour(x,y,psi,[-.4 -.4],'g','LineWidth',1.1);

p4=contour(x,y,psi,[-.6 -.6],'b','LineWidth',1.1);

p5=contour(x,y,psi,[-.8 -.8],'c','LineWidth',1.1);

p6=contour(x,y,psi,[-1 -1],'m','LineWidth',1.1);

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$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

%ylabel({'$\eta=1\quad$'},'Interpreter','latex','FontSize',10,'rot',360,'FontName','Times New Roman','FontWeight','Normal');

%title('Happel$^\prime$s model','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

axis square

axis on

[Shreen El-Sapa]

Why the program excuted the final, p6 only?

Answer 1

[VBBV]

Use hold on as shown below

beta1=0.00001; beta2=0.00001;
Zeta1=1;Zeta2=1;
k=.00001;U=1;
alpha1=real(sqrt(Zeta1.^2./2+Zeta1.*sqrt(Zeta1.^2-4.*k.^2)./2));
alpha2=real(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 = -(Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 4 .* alpha2 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 4 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .* beta1 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .* alpha2 .^ 2 .* beta1 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 3 .* alpha2 .^ 2 .* beta1 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 2 .* alpha2 .^ 3 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 2 .* beta2 + 2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 4 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 5 .* alpha1 .^ 4 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .* beta1 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 .* beta2 + 4 .* Zeta1 .^ 4 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 2 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 3 + 20 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 2 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .* alpha2 .^ 3 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 4 + 8 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 4 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 4 .* alpha1 .^ 3 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta2 + 12 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta2 + 3 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 .* beta2 + 8 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 12 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 4 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* beta1 + 40 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 36 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 + 36 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 3 .* beta1 + 4 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 4 .* alpha2 .^ 2 + 8 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 4 + 8 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 4 .* alpha2 .^ 2 + 16 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 4 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .* alpha1 .* beta2 + 6 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .* alpha2 .* beta2 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 12 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .^ 2 .* beta1 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 24 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .^ 2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 + 72 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta1 + 72 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta1 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta1 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 .^ 2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .^ 3 + 24 .* Zeta1 .^ 4 .* a .* alpha1 .^ 2 + 48 .* Zeta1 .^ 4 .* a .* alpha1 .* alpha2 + 12 .* Zeta1 .^ 4 .* a .* alpha1 .* beta1 + 24 .* Zeta1 .^ 4 .* a .* alpha2 .^ 2 + 12 .* Zeta1 .^ 4 .* a .* alpha2 .* beta1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 + 96 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* beta1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 + 24 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .* beta1 + 24 .* Zeta1 .^ 4 .* alpha1 + 24 .* Zeta1 .^ 4 .* alpha2 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* alpha1 + 48 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* alpha2) .* a ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* beta2 + 2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 2 + 4 .* Zeta1 .^ 4 .* a .^ 2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 + 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 + 8 .* Zeta1 .^ 4 .* a + 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 + 8 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 + 4 .* Zeta1 .^ 2 .* a .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 + 16 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 + 8 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 + 4 .* Zeta1 .^ 2 .* alpha1 + 4 .* Zeta1 .^ 2 .* alpha2 + 8 .* Zeta2 .^ 2 .* alpha1 + 8 .* Zeta2 .^ 2 .* alpha2) ./ alpha2 .^ 2 ./ alpha1 .^ 2;
BB = -0.3e1 .* exp(alpha1 .* a) .* sqrt(alpha1 .* a) .* (Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* beta2 + 0.4e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha2 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha2 + 0.4e1 .* Zeta1 .^ 2 + 0.8e1 .* Zeta2 .^ 2) .* sqrt(a) .* Zeta1 .^ 2 .* (a .* beta1 + 0.2e1) .* sqrt(0.2e1) ./ alpha1 .* pi .^ (-0.1e1 ./ 0.2e1) ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 .* beta2 - Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 .* beta2 + 0.2e1 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 - 0.2e1 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 + 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 - 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 .* beta2 + 0.4e1 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 - 0.4e1 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 + 0.4e1 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* beta1 - 0.4e1 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .* beta1 + 0.12e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 0.8e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 - 0.12e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 - 0.8e1 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 .* beta2 + 0.4e1 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 0.4e1 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 0.8e1 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 0.8e1 .* Zeta1 .^ 4 .* a .* alpha1 - 0.8e1 .* Zeta1 .^ 4 .* a .* alpha2 + 0.16e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 - 0.16e2 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .^ 3 + 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 0.4e1 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* a .* alpha2 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .^ 3 + 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 0.8e1 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 0.8e1 .* Zeta2 .^ 2 .* a .* alpha2 .^ 3 + 0.4e1 .* Zeta1 .^ 2 .* alpha1 .^ 2 - 0.4e1 .* Zeta1 .^ 2 .* alpha2 .^ 2 + 0.8e1 .* Zeta2 .^ 2 .* alpha1 .^ 2 - 0.8e1 .* Zeta2 .^ 2 .* alpha2 .^ 2);
CC =0.3e1 .* (Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 4 .* Zeta1 .^ 2 .* a .* alpha1 + 8 .* Zeta2 .^ 2 .* a .* alpha1 + 4 .* Zeta1 .^ 2 + 8 .* Zeta2 .^ 2) .* sqrt(a) .* sqrt((alpha2 .* a)) .* (Zeta1 .^ 2) .* (a .* beta1 + 2) .* exp((alpha2 .* a)) .* sqrt(0.2e1) ./ alpha2 .* pi .^ (-0.1e1 ./ 0.2e1) ./ (Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* beta1 .* beta2 - Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .* beta1 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 3 .* beta2 + 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta2 - 2 .* Zeta1 .^ 4 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 + 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* beta1 - 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .* alpha2 .^ 2 - 4 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 3 .* alpha2 .^ 2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* beta2 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 .* alpha2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 3 .* beta2 + 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 3 .* alpha2 .^ 2 - 4 .* Zeta2 .^ 2 .* a .^ 3 .* alpha1 .^ 2 .* alpha2 .^ 3 + 4 .* Zeta1 .^ 4 .* a .^ 2 .* alpha1 .* beta1 - 4 .* Zeta1 .^ 4 .* a .^ 2 .* alpha2 .* beta1 + 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 2 + 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* beta1 - 12 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .^ 2 - 8 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .^ 2 .* alpha2 .* beta1 + Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* beta2 - Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 .^ 2 .* beta2 + 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 4 .* Zeta1 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .^ 3 .* alpha2 - 8 .* Zeta2 .^ 2 .* a .^ 2 .* alpha1 .* alpha2 .^ 3 + 8 .* Zeta1 .^ 4 .* a .* alpha1 - 8 .* Zeta1 .^ 4 .* a .* alpha2 + 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha1 - 16 .* Zeta1 .^ 2 .* Zeta2 .^ 2 .* a .* alpha2 + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 3 + 4 .* Zeta1 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 4 .* Zeta1 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 4 .* Zeta1 .^ 2 .* a .* alpha2 .^ 3 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 3 + 8 .* Zeta2 .^ 2 .* a .* alpha1 .^ 2 .* alpha2 - 8 .* Zeta2 .^ 2 .* a .* alpha1 .* alpha2 .^ 2 - 8 .* Zeta2 .^ 2 .* a .* alpha2 .^ 3 + 4 .* Zeta1 .^ 2 .* alpha1 .^ 2 - 4 .* Zeta1 .^ 2 .* alpha2 .^ 2 + 8 .* Zeta2 .^ 2 .* alpha1 .^ 2 - 8 .* Zeta2 .^ 2 .* alpha2 .^ 2);
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

psi = 201×201

-1.0699 -1.0740 -1.0781 -1.0822 -1.0864 -1.0905 -1.0946 -1.0988 -1.1029 -1.1071 -1.1112 -1.1154 -1.1196 -1.1238 -1.1279 -1.1321 -1.1363 -1.1405 -1.1446 -1.1488 -1.1530 -1.1572 -1.1613 -1.1655 -1.1696 -1.1738 -1.1780 -1.1821 -1.1862 -1.1903 -1.0526 -1.0567 -1.0608 -1.0649 -1.0690 -1.0731 -1.0772 -1.0813 -1.0854 -1.0895 -1.0937 -1.0978 -1.1020 -1.1061 -1.1102 -1.1144 -1.1186 -1.1227 -1.1269 -1.1310 -1.1352 -1.1393 -1.1435 -1.1476 -1.1517 -1.1559 -1.1600 -1.1641 -1.1683 -1.1724 -1.0354 -1.0395 -1.0435 -1.0476 -1.0516 -1.0557 -1.0598 -1.0639 -1.0680 -1.0721 -1.0762 -1.0803 -1.0844 -1.0885 -1.0926 -1.0967 -1.1009 -1.1050 -1.1091 -1.1133 -1.1174 -1.1215 -1.1256 -1.1298 -1.1339 -1.1380 -1.1421 -1.1462 -1.1503 -1.1544 -1.0183 -1.0223 -1.0263 -1.0303 -1.0343 -1.0384 -1.0424 -1.0465 -1.0505 -1.0546 -1.0587 -1.0628 -1.0669 -1.0710 -1.0750 -1.0791 -1.0833 -1.0874 -1.0915 -1.0956 -1.0997 -1.1038 -1.1079 -1.1120 -1.1161 -1.1202 -1.1243 -1.1284 -1.1324 -1.1365 -1.0012 -1.0052 -1.0091 -1.0131 -1.0171 -1.0211 -1.0251 -1.0292 -1.0332 -1.0372 -1.0413 -1.0453 -1.0494 -1.0535 -1.0575 -1.0616 -1.0657 -1.0698 -1.0738 -1.0779 -1.0820 -1.0861 -1.0902 -1.0943 -1.0983 -1.1024 -1.1065 -1.1105 -1.1146 -1.1187 -0.9842 -0.9881 -0.9921 -0.9960 -1.0000 -1.0039 -1.0079 -1.0119 -1.0159 -1.0199 -1.0239 -1.0280 -1.0320 -1.0360 -1.0401 -1.0441 -1.0482 -1.0522 -1.0563 -1.0603 -1.0644 -1.0685 -1.0725 -1.0766 -1.0806 -1.0847 -1.0887 -1.0928 -1.0968 -1.1009 -0.9672 -0.9711 -0.9750 -0.9789 -0.9829 -0.9868 -0.9908 -0.9947 -0.9987 -1.0027 -1.0067 -1.0107 -1.0147 -1.0187 -1.0227 -1.0267 -1.0307 -1.0347 -1.0388 -1.0428 -1.0468 -1.0509 -1.0549 -1.0589 -1.0630 -1.0670 -1.0710 -1.0751 -1.0791 -1.0831 -0.9503 -0.9542 -0.9581 -0.9620 -0.9659 -0.9698 -0.9737 -0.9776 -0.9815 -0.9855 -0.9894 -0.9934 -0.9974 -1.0013 -1.0053 -1.0093 -1.0133 -1.0173 -1.0213 -1.0253 -1.0293 -1.0333 -1.0374 -1.0414 -1.0454 -1.0494 -1.0534 -1.0574 -1.0614 -1.0654 -0.9335 -0.9373 -0.9412 -0.9450 -0.9489 -0.9528 -0.9567 -0.9605 -0.9645 -0.9684 -0.9723 -0.9762 -0.9802 -0.9841 -0.9881 -0.9920 -0.9960 -0.9999 -1.0039 -1.0079 -1.0119 -1.0159 -1.0199 -1.0238 -1.0278 -1.0318 -1.0358 -1.0398 -1.0438 -1.0478 -0.9168 -0.9206 -0.9244 -0.9282 -0.9320 -0.9359 -0.9397 -0.9436 -0.9474 -0.9513 -0.9552 -0.9591 -0.9630 -0.9669 -0.9708 -0.9748 -0.9787 -0.9826 -0.9866 -0.9905 -0.9945 -0.9985 -1.0024 -1.0064 -1.0103 -1.0143 -1.0183 -1.0222 -1.0262 -1.0302

hold on
%[DH,h2]=contour(x,y,psi,5,'-k','ShowText','on');
p1=contour(x,y,psi,[-.1 .1],'k','LineWidth',1.1);  %,'ShowText','on'
p2=contour(x,y,psi,[-.2 .2],'r','LineWidth',1.1);  
p3=contour(x,y,psi,[-.4 .4],'g','LineWidth',1.1);  
p4=contour(x,y,psi,[-.6 .6],'b','LineWidth',1.1);  
p5=contour(x,y,psi,[-.8 .8],'c','LineWidth',1.1);  
p6=contour(x,y,psi,[-1 1],'m','LineWidth',1.1);  

m1=100;
r1=ones(1,m1+1)*a;
th=[0:2*pi/m1:2*pi];
phi = 0:pi/50:2*pi;
polar(th,r1,'k');
fill(a*cos(phi),a*sin(phi),'y');
title('$\kappa=0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')
%ylabel({'$\eta=1\quad$'},'Interpreter','latex','FontSize',10,'rot',360,'FontName','Times New Roman','FontWeight','Normal');
%title('Happel$^\prime$s model','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')
axis square
axis on