beta11 = @(z) 1/2*(1+3*z/2+(1+z+9/4*z.^2).^(0.5)); beta12 = @(z) 1/2*(1+3*z/2-(1+z+9/4*z.^2).^(0.5)); beta21 = @(z) 1/2./(1-5/12*z).*(1+2*z/3+(1+z+7/12*z.^2).^(0.5)); beta22 = @(z) 1/2./(1-5/12*z).*(1+2*z/3-(1+z+7/12*z.^2).^(0.5)); [x,y] = meshgrid([-4:0.005:4]); z = x+i*y; p1 = abs(beta11(z)); p2 = abs(beta12(z)); p3 = abs(beta21(z)); p4 = abs(beta22(z)); subplot(2,1,1); imagesc([-4,4],[-4,4],(p1<1).*(p2<1)) axis equal xlim([-4 4]); title('A.B. 2-step') subplot(2,1,2); imagesc([-4,4],[-4,4],(p3<1).*(p4<1)) axis equal xlim([-4 4]); title('A.M. 2-step')