%% 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