%% Plot of A-sability regions
% setting theta values, gridspace and grid nodes
theta = [0,0.2,0.4,0.48,0.5,0.52,0.6,0.8,1];
dx = 0.01;
[x,y] = meshgrid(-4:dx:4);
figure()
% loop for plotting stability domain for each theta
for ii =1:length(theta)
% If |R(z)|<1 then 1 else 0
R = (1-heaviside(2*x+(1-2*theta(ii))*(x.^2+y.^2)));
% subplotting the domain
subplot(3,3,ii);
p(ii) = pcolor(x,y,R); shading flat;
xlabel('Re(z)')
ylabel('Im(z)')
title(['Red marked stability domain for theta =',num2str(theta(ii))])
end