>speryr=31556926;
>cmperpc=3.086e18;
>alphab=2.59e-13 * speryr/(cmperpc**3);
>q=5e48*speryr;
>n=10*(cmperpc**3);
>function drdt(t,r)
$useglobal;
$return q/(4*3.14159265*r*r*n) - alphab*n*r/3;
$endfunction
>function norecombr(t)
$useglobal;
$return (3*q*t/(4*3.14159265*n))**(1/3);
$endfunction
>t=0:1:6e4;
>samp=1:6e3:6e4;
>xplot(t,norecombr(t));
>hold on
>r=runge("drdt",t,0.1);
>r[samp]
>xplot(t,r);