Simple modelling of seasonal variation of temperature? I'm really curious about this:
What is the simplest (or most simplified) differential equation that accounts for the variations of temperature throughout the year at some point on the northern hemisphere?
It should explain why, for example, even though in the weeks leading up to the Fall equinox, when days are still longer than nights, temperature is slowly getting cooler.
I don't know anything about physics, I work in mathematics, so please forgive me if this question is too ambiguous. 
 A: I'd look at it as an energy storage vs loss situation.
Take a patch of earth (square slab) and neglect rotation of the earth around its axis (days) so that the patch always faces the sun.  At any time it's receiving an incident solar flux (assume constant) and emitting due to its own temperature.  The slab also has some thermal mass (capturing the ground, air, water etc).
Where the seasons come about due to the orientation of the earth relative to the sun. Because axis of rotation of the earth  is not parallel to that of the orbit, your patch of earth will vary between being pointed fully at the sun to being angled a bit away. That angling reduces the cross section of the patch of earth relative the to solar flux, reducing the total energy absorbed (imagine the extreme case of turning it 90 deg st all of the solar rays are parallel to the surface).
Now let's get a little more mathy.  Let's ignore all the real geometry of the situation and just say that the cross section of the patch relative to the incident flux varies as a simple $cos$ st.
$$q''_{incident}=q''_0(1-C(1+cos(2\pi\:t))$$
I think that we can get away with calling the loss term  a linear with $q''_l$.  It should vary as $T^4$, but meh.
$$ q''_l = q''_{l,0}T$$
Last we put in a storage term to track the rate of temperature change. It's mass per unit area of our slab $m''$, specific heat $C_p$ and $\dot T$.  Put them all together and get
$$m''C_p \dot T = q''_0(1-C(1+cos(2\pi\:t)) - q''_{l,0}T$$
here are numerical results with $m''C_p=30$, $C=0.2$, $q''_{i,0}=q''_{l,0}=2$. The lag that you see between incident and loss is due to the thermal inertia of the system. (apologies for the excel plot)

