q1 = @(z) 1+z;
q2 = @(z) 1+z+z.^2/2;
q3 = @(z) 1+z+z.^2/2+z.^3/6;
q4 = @(z) 1+z+z.^2/2+z.^3/6+z.^4/24;
[x,y] = meshgrid([-4:0.005:4]);
z = x+i*y;
p1 = abs(q1(z));
p2 = abs(q2(z));
p3 = abs(q3(z));
p4 = abs(q4(z));
subplot(2,2,1);
imagesc([-4,4],[-4,4],p1<1)
axis equal
xlim([-4 4]);
title('Euler')

subplot(2,2,2);
imagesc([-4,4],[-4,4],p2<1)
axis equal
xlim([-4 4]);
title('2nd Order')

subplot(2,2,3);
imagesc([-4,4],[-4,4],p3<1)
axis equal
xlim([-4 4]);
title('3rd Order')

subplot(2,2,4);
imagesc([-4,4],[-4,4],p4<1)
axis equal
xlim([-4 4]);
title('4th Order')