This experiment (with tea and milk) seems to show that immediate addition of the milk gives the warmest drink:
But the experiment is far from perfect, especially because of the lack of replication.
And this simplified derivation, based on Newton's Law of Cooling finds the same. But it too is very open to well-reasoned critiques. It certainly relies on one rather questionable assumption.
And here's my own experiment. Experimental protocol is summarised at the end of this post.
I tracked both methods of drinking coffee over a $15\mathrm{min}$ period of time. At the end I found the coffee with immediate addition of milk to be about $8^\circ$$\mathrm{Celsius}$ hotter than when the milk was added at the end of the run. Even without replication that's likely to be highly significant.
The dimensionless temperature $\Theta$ v. time plots are quite revealing, Series $1$ is with immediate addition of milk, Series $2$ with addition at the end:
Assuming Newton's Law of cooling (see link for symbols used) applies, we can write:
$$\ln \Theta=\ln\Big[\frac{T(t)-T_{\infty}}{T_0-T_{\infty}}\Big]=-\frac{t}{\tau}$$
where:
$$\frac{1}{\tau}=\frac{h A}{m c_p}$$
$\tau$ is the so-called characteristic time. For immediate addition of the milk, Linear Regression gave a value of $1/\tau=0.0123\text{ 1/s}$ and $1/\tau=0.0148\text{ 1/s}$ for addition at the end, from the two runs above. That's a difference of about $17\text{ %}$.
Firstly, strong linearity of the $(\ln\Theta, t)$ plots indicates that radiative and evaporative losses are probably negligible compared to convective losses. Newton's Law of cooling seems to hold.
Secondly, the greater cooling rate for Series $2$ indicates a higher value for $\frac{1}{\tau}$, so for $\frac{hA}{mc_p}$.
A possible cause of the difference in cooling rate may lay in the factor $A/m$. For a cylinder with (constant) radius $R_0$ we found:
$$\frac{A}{m}=2\pi \Big(\frac{R_0^2}{m}+\frac{1}{\pi \rho R_0}\Big)$$
As $m$ goes up, $A/m$ comes down and so does $1/\tau$ and the cooling rate. Surely a bucket of hot coffee would take longer to cool than a mere cup?
So maybe what we're observing is a mere mass effect.
The influence of the $A/m$ ratio on cooling rate was investigated empirically as follows.
The first two experiments were replicated but using coffee only. Another run was added in between the first two so we have three cooling rates for $A/m=0.83$, $A/m=0.87$, $A/m=0.91$:
Here the effect of $A/m$, although still present, is much smaller than in the first run.
Could the larger difference in cooling rate be caused by using milk? Milk has a lower heat capacity than water because it contains fat. That would increase $1/\tau$.
Various internet sources put the $c_p$ value for semi-skimmed milk at about $4.0\text{ kJ/(kgK)}$, as opposed to $4.18\text{ kJ/(kgK)}$ for pure water, a difference of about $4\text{ %}$. But we only used $27\text{ %}$ of milk in our first two runs, so it's unlikely to account for the $17\text{ %}$ observed difference in $1/\tau$ values.
Experiment protocol (summary)
$230\mathrm{g}$ of black, boiling instant coffee with $27\text{ %}$ of semi-skimmed milk.
Paper cup (take away), no lid. No insulation.
Slow magnetic stirring.
Environmental temperature $T_{\infty}=19.0^\circ\mathrm{Celsius}$, for the first two runs. For later runs that value had to be adjusted slightly.
