I wondered about this myself and decided to sit down and do it today out of curiosity. Post is 2 years old but it's worth putting. Also it's my first time using Latex so forgive the inconsistency of lettering.
Start from Newton's law of cooling, with t as time, T temp, and k constant,
\frac{dT_{a}}{dt} = -ka(T_{a}-T_{b})
Which we picked the sign of to make sure if T_{a} >T_{b}, that T will decrease.
Now, we remember the definition of specific heat, C_{a}=\frac{dU_{a}}{dT}, then
\frac{dU_{a}}{dt}=\frac{dU_{a}}{dT}*\frac{dT_{a}}{dt}=C*\frac{dT_{a}}{dt},
So our new Newton's law is \frac{dU_{a}}{dt} = -\frac{k_{a}}{C_{a}}(T_{a}-T_{b})= -k(T_{a}-T_{b})
And now we have to choose a T dependence for U. For Newton's law, we implicitly assumed U proportional to T, as is true for both the ideal gas and for metals from Dulong-Petit rule, so then we can express U=C*T, and separate variables for integration;
∫\frac{dU_{a}}{\frac{U_{a}}{C_{a}}-\frac{U_{b}}{C_{b}}}=∫-k*dt
\frac{Ln[{\frac{U_{a}}{C_{a}}-\frac{U_{b}}{C_{b}}}]}{\frac{1}{C_{a}}-\frac{\frac{dU_{b}}{dU_{a}}}{C_{b}}}=-k*t+c,
c integration constant for later. Now we employ energy conservation, so that
U_{a}+U_{b}=E, E constant. therefore 1+\frac{dU_{b}}{dU_{a}}=0, so our equation becomes;
Ln[{\frac{U_{a}}{C_{a}}-\frac{U_{b}}{C_{b}}}]=(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t+c)=(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t)+d
Exponentiate, and absorb the d as multiplicative integration constant. Add over U_{b}, multiply by C_{a}, and you get;
U_{a}=C_{a}*A*e^{(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t)}+\frac{C_{a}}{C_{b}}U_{b}
And to separate the coupled equations, we again use conservation of energy, so
U_{b}=E-U_{a}, add over the U term, factor out U, divide by coefficient, and;
U_{a}=\frac{C_{a}*A*e^{(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t)}+\frac{C_{a}}{C_{b}}E}{1+\frac{C_{a}}{C_{b}}}
And sub in U_{a}=C_{a}T_{a}, so divide whole equation by Ca,
T_{a}=\frac{A*e^{(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t)}+\frac{E}{C_{b}}}{1+\frac{C_{a}}{C_{b}}},
With E=C_{a}T_{a}+C_{b}T_{b}=C_{a}T_{a0}+C_{b}T_{b0}, since E same at all times.
Then we apply initial conditions, at t=0 T_{a}=T_{a0}=\frac{A+\frac{C_{a}T_{a0}+C_{b}T_{b0}}{C_{b}}}{1+\frac{C_{a}}{C_{b}}},
So A=T_{a0}-T_{b0} and we have our final form;
T_{a}=\frac{(T_{a0}-T_{b0})e^{(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t)}+\frac{C_{a}T_{a}+C_{b}T_{b}}{C_{b}}}{1+\frac{C_{a}}{C_{b}}}
And we check some limits for sanity; at t→∞,
T_{a}=\frac{C_{a}T_{a0}+C_{b}T_{b0}}{C_{a}+C_{b}},
So if C_{b} goes to 0, ie you get no heat from changing it's temperature, then T_{a}=T_{a0} ; it doesn't change temperature since there is no energy in the other system
If C_{b} goes to ∞, it takes infinite energy to change it's temperature, and we see by L'Hopitals rule that T_{a}=T_{b0}; it raises to the unmovable temperature
And finally if both C's are equal, T_{a}=\frac{T_{a0}+T_{b0}}{2}; it reaches equilibrium at the average of the initial temperatures.
Nonlinear U dependences on T will need to be solved differently from the step listed as such, but for first order and most common systems this should work. Hope it helps.
