How I can perform SIMO channel convolution in MATLAB
Your question is specific to how to do the convolution with a matrix multiplication in the case of a SIMO system. 

Below is a graphic from Compensating Loudspeaker frequency response in an audio signal showing how the Toeplitz matrix is used to do SISO convolution with matrix multiplication:



In your case "r" in the graphic is "x" in your code, and the Toeplitz Matrix A is formed by placing x in the column and shifting this along the diagonal over the number of columns given by your channel length (in your case 12). So if your transmit signal was 160 samples long and your channel 16 long, A would be a 175 x 16 matrix. As shown in the figure above, if you only had a single channel, the channel would be multiplied as a column vector resulting in the matrix product of 175x16 with 16x1 resulting in a solution as a column vector of 175 x 1. (Remember to do the matrix multiplication, the inside dimensions must match and the result will be the outside dimensions).

For multiple channels as your case, simply multiply by the 16x4 channel matrix instead. This will result in 175x16 mulitplied with 16x4 resulting in a matrix result of 175 x 4, representing a column for each received signal.