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

