function dYds = deriv(x,y)
global t kappa alpha V h rho r
d2Tds2 = ( h*y(1)/t - V^2/(4*pi^2*r^2*rho)/(1+alpha*y(1)) ) / kappa;
dYds = [y(2); d2Tds2];