Tv=1000; % K
T0=300; % K
Q=1.38e-23;
L=1e-4; % m.
syms T(t)
DE = diff(T) == -Q*(Tv^4-T0^4)/3*L*T;
DEs = dsolve(DE, T(0)==T0)
figure
fplot(DEs, [0 1E+16])
grid
xlabel('$t$', 'Interpreter','latex')
ylabel('$T(t)$', 'Interpreter','latex')
title(['$T(t) = ' latex(DEs) '$'], 'Interpreter','latex')
