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,
[itex]\frac{dT_{a}}{dt}[/itex] = -ka(T[itex]_{a}[/itex]-T[itex]_{b}[/itex])
Which we picked the sign of to make sure if T[itex]_{a}[/itex] >T[itex]_{b}[/itex], that T will decrease.
Now, we remember the definition of specific heat, C[itex]_{a}[/itex]=[itex]\frac{dU_{a}}{dT}[/itex], then
[itex]\frac{dU_{a}}{dt}[/itex]=[itex]\frac{dU_{a}}{dT}[/itex]*[itex]\frac{dT_{a}}{dt}[/itex]=C*[itex]\frac{dT_{a}}{dt}[/itex],
So our new Newton's law is [itex]\frac{dU_{a}}{dt}[/itex] = -[itex]\frac{k_{a}}{C_{a}}[/itex](T[itex]_{a}[/itex]-T[itex]_{b}[/itex])= -k[itex](T_{a}-T_{b})[/itex]
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;
∫[itex]\frac{dU_{a}}{\frac{U_{a}}{C_{a}}-\frac{U_{b}}{C_{b}}}[/itex]=∫-k*dt
[itex]\frac{Ln[{\frac{U_{a}}{C_{a}}-\frac{U_{b}}{C_{b}}}]}{\frac{1}{C_{a}}-\frac{\frac{dU_{b}}{dU_{a}}}{C_{b}}}[/itex]=-k*t+c,
c integration constant for later. Now we employ energy conservation, so that
[itex]U_{a}+U_{b}[/itex]=E, E constant. therefore 1+[itex]\frac{dU_{b}}{dU_{a}}[/itex]=0, so our equation becomes;
[itex]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[/itex]
Exponentiate, and absorb the d as multiplicative integration constant. Add over [itex]U_{b}[/itex], multiply by [itex]C_{a}[/itex], and you get;
[itex]U_{a}=C_{a}*A*e^{(\frac{1}{C_{a}}+\frac{1}{C_{b}})(-k*t)}+\frac{C_{a}}{C_{b}}U_{b}[/itex]
And to separate the coupled equations, we again use conservation of energy, so
[itex]U_{b}=E-U_{a}[/itex], add over the U term, factor out U, divide by coefficient, and;
[itex]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}}}[/itex]
And sub in [itex]U_{a}=C_{a}T_{a}[/itex], so divide whole equation by Ca,
[itex]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}}}[/itex],
With [itex]E=C_{a}T_{a}+C_{b}T_{b}=C_{a}T_{a0}+C_{b}T_{b0}[/itex], since E same at all times.
Then we apply initial conditions, at t=0 [itex]T_{a}=T_{a0}=\frac{A+\frac{C_{a}T_{a0}+C_{b}T_{b0}}{C_{b}}}{1+\frac{C_{a}}{C_{b}}}[/itex],
So [itex]A=T_{a0}-T_{b0}[/itex] and we have our final form;
[itex]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}}}[/itex]
And we check some limits for sanity; at t→∞,
[itex]T_{a}=\frac{C_{a}T_{a0}+C_{b}T_{b0}}{C_{a}+C_{b}}[/itex],
So if [itex]C_{b}[/itex] goes to 0, ie you get no heat from changing it's temperature, then [itex]T_{a}=T_{a0}[/itex] ; it doesn't change temperature since there is no energy in the other system
If [itex]C_{b}[/itex] goes to ∞, it takes infinite energy to change it's temperature, and we see by l'hospital's rule that [itex]T_{a}=T_{b0}[/itex]; it raises to the unmovable temperature
And finally if both C's are equal, [itex]T_{a}=\frac{T_{a0}+T_{b0}}{2}[/itex]; 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.