Do you want a single temperature for each decade, or you you want an average for each lat/lon coordinate pair?
I wasn't able to get splitapply or groupsummary to work with a 3D dataset, but I was able to achieve group means with a for-loop.
The challenge here is perhaps identifying the groups. I want to convert your time variable to a datetime so I can group by year. The following code allows we to inspect the attributes of time.
info = ncinfo('keeganCarvalho_bob.nc')
From here, I can see that time is the number of days since 1800-1-1 00:00:00, and that 
is 1 day. I performed the convertion with the following code.
is 1 day. I performed the convertion with the following code.dTime = datetime(1800,1,1)+days(time)
I then extracted the year and used the discretize function to group the years into your specified bins.
yr = year(dTime);
G = discretize(yr,[1986 1996 2006 2016 2019]);
All that is left is to get the mean.
for d = 1:max(G)
decadeMean(d) = mean(temp(:,:,G==d),"all");
end
plot(decadeMean)
If instead you want the mean for each lat/lon location, you could try something like this. You'd have to change the plotting code to be something appropriate for matrix data.
for d = 1:max(G)
decadeMean(:,:,d) = mean(temp(:,:,G==d),3);
end
% rearrange data so decade means are in rows, and columns are lat/lon pairs
plot(reshape(permute(decadeMeanLL,[3,1,2]),max(G),[]))
