I don't see anything wrong with Akira's approach. It does what is asked (except that it has a wrong start of season).
However, I'd use a different approach. I'd simply reshape into an Mx3xS (where M is data points per month and S is number of seasons) and sum across the columns. You may need to add some extra months full of 0 (which don't affect the sum) to make sure that we've got 3 months for all seasons.
So, if the matrix starts in January and ends in December, we need to add one month of 0s at the start and 2 months of 0s at the end:
padded = cat(3, zeros(size(data, 1), size(data, 2), 'like', data), data);
padded = cat(3, padded, zeros(size(data, 1), size(data, 2), mod(-size(padded, 3), 3), 'like', data));
seasonal = reshape(padded, size(data, 1) * size(data, 2), 3, );
seasonalsum = sum(seasonal, 2);
seasonalsum = reshape(seasonalsum, size(data, 1), size(data, 2), );
edit: Note that I'm padding with 0s since it doesn't affect the sum. If you were to calculate the mean or some other statistics, I'd pad with NaNs instead, and use the 'omitnan' option of the corresponding function. You could do that for the sum as well.