Step 6:  

Add a pooling layer

clear all;
close all;
clc;

[XTrain,YTrain] = digitTrain4DArrayData;
size(XTrain)      %images
size(YTrain)      %correct answer labels

XTrain=1-XTrain;  % Reverse the black and white colors.  Save and run the program to see the difference.

perm = randperm(size(XTrain,4),20);  % Randomize the order of images in XTrain
for i = 1:20
subplot(4,5,i);
imshow(XTrain(:,:,:,perm(i)));
end
 
layers = [
imageInputLayer([28 28 1])   
convolution2dLayer(3,8,'Padding','same')
reluLayer
maxPooling2dLayer(2,'Stride',2)

We'll just go with the default values of the reluLayer (this just zeros out all the negative numbers in the filtered image, and leaves the rest unchanged).
For the pooling layer, we're deciding to use a 2x2 compression (so every 2x2 patch of the original image will get reduced to one pixel wihich is the maximum value of the four pixels).  We also will use a stride of 2, which means that after we pool one set of 2x2 pixels, we'll jump over by 2 pixels to the right and repeat the process).