Оптимизация кода

clear; clc; close;
 
d = 0.1; ai = -0.95:0.01:3.15; bi = -1.05:0.01:1.25; [a,b] = meshgrid(ai,bi);
 
l1 = b < d*a+1-d;
l2 = b > (a-2)*(1+d);
l3 = b > (d-1)*a;
l = l1&l2&l3;
C = contour(a,b,l,1,'r');
 
gn = 0.5;
num_ab = 0;
num_fi = 0;
num_x = 0;
HeycT = true;
kx = -2:4/50:2;
kfi = -0.99*pi:1.98*pi/50:0.99*pi;
maccuB_fi = zeros(length(kx)*length(kfi),1);
maccuB_x = maccuB_fi;
maccuB_ab = zeros(length(ai)*length(bi),2);
R = eps;
tic
for a = ai
    a
    for b = bi
        if (b < d*a+1-d) && (b > (a-2)*(1+d)) && (b > (d-1)*a)          
            for x1 = kx
                for fi1 = kfi
                    for n = 1:2000
                        fi2 = fi1-a*fi1+x1;
                        x2 = -b*fi1+d*x1+gn;
                        while or(fi2 > pi, fi2 < -pi) == 1
                            if fi2 > pi
                                fi2 = fi2-2*pi;
                            elseif fi2 < -pi
                                fi2 = fi2+2*pi;
                            end;
                        end;
                        fi1 = fi2;
                        x1 = x2;
                        if a-d*a+b == 0
                            x_cT = a*gn/eps;
                            fi_cT = gn/eps;
                        else
                            x_cT = a*gn/(a-d*a+b);
                            fi_cT = gn/(a-d*a+b);
                        end;
                        if (fi2-fi_cT)^2/R^2 + (x2-x_cT)^2/R^2 < 1
                            HeycT = false;
                            break;
                        elseif n == 2000
                            HeycT = true;
                        end;
                    end;
                    num_fi = num_fi+1;
                    num_x = num_x+1;
                    maccuB_fi(num_fi,1) = fi2;
                    maccuB_x(num_x,1) = x2;
                    if HeycT == true
                        break;
                    end;         
                end;
                if HeycT == true
                    break;
                end;
            end;
            if HeycT == false
                num_ab = num_ab+1;
                maccuB_ab(num_ab,1) = a;
                maccuB_ab(num_ab,2) = b;
            end;
            HeycT = true;
            num_fi = 0;
            num_x = 0;
            maccuB_fi = zeros(length(kx)*length(kfi),1);
            maccuB_x = maccuB_fi;
        end;
    end;
end;
toc
hold on;
grid on;
plot(maccuB_ab(:,1), maccuB_ab(:,2), '.b');
 
xlabel('a'); ylabel('b');

i = 0;
for a = ai
    for b = bi
        if (b < d*a+1-d) && (b > (a-2)*(1+d)) && (b > (d-1)*a)
            i = i+1;
            BeKTop_a(i,1) = a;
            BeKTop_b(i,1) = b;
        end;
    end;
end;
 
maccuB_ab = zeros(length(BeKTop_a),2);
 
a = BeKTop_a;
b = BeKTop_b;