Let's start with a row vector of January 1sts.
Jan1st = datetime(2010:2021, 1, 1);
We want to compute a vector of days for each January 1st that increases by 8 days until the end of the year. Luckily the number of days in the year is not a multiple of 8, so each year will have the same number of values represented.
dayIncrements = days(0:8:365).';
If you were using release R2020b or later you could use implicit expansion to sum Jan1st and dayIncrements to generate a matrix. If you need this as a vector you could reshape it. Since you indicated you're using release R2018a you'd need to repmat the two vectors.
t = Jan1st + dayIncrements;
t2 = repmat(Jan1st, numel(dayIncrements), 1) + repmat(dayIncrements, 1, numel(Jan1st));
check = isequal(t, t2)
Now let's look at the first and last days in each column. Leap years end on December 26th, non leap years on December 27th.
firstAndLast = t([1 end], :)
