Ошибка в дифференциальном уравнении

//Параметры нагревателей в количестве 10 шт
U=380;//Напряжение сети переменного трехфазного тока, В
nn=10;//Количество нагревателей
mn=10;//Масса нагревателя, кг
dn=0.06;//Диаметр нагревателя, м
ln=1;//Длина нагревателя, м
Cn=1;//Удельная теплоемкость нагревателя, кДж/(кг*К)
En=0.65;//Степень черноты нагревателя
//Габаритные размеры печи, м
A=0.6;//Ширина
B=1.1;//Длина
H=0.6;//Высота
bb=0.1;//Толщина стенки
//Размеры окна, м
a=0.3;//Ширина
b=0.2;//Высота
//Физические параметры стен печи
rost=2100;//Плотность материала стенки печи, кг/м^3
Cst=1;//Удельная теплоемкость стенки печи, кДж/(кг*К)
lst=1.5;//Коэффициент теплопроводности, Вт/(м*К)
//Параметры заготовок в количестве 5 шт
nz=10;//Количество заготовок
mz=12;//Масса заготовки, кг
Cz=0.8;//Удельная теплоемкость заготовки, кДж/(кг*К)
//Дополнительные величины
ast=50;//Коэффициент теплоотдачи от наружных поверхностей печи в о.с., Вт/(м^2*К)
tst2=45;//Температура наружной поверхности печи, оС
Eok=0.8;//Степень черноты окна
tos=20;//Температура окружающей среды, оС
//Расчеты
//Площадь окна, м^2
Fok=a*b;
//Площадь внутренних поверхностей печи, м^2
F1=2*((B-2*bb)*(A-2*bb)+(H-2*bb)*(A-2*bb)+(B-2*bb)*(H-2*bb)-a*b);
//Площадь наружных поверхностей печи, м^2
Fnar=2*(B*A+H*A+B*H-a*b);
//Площадь поверхности нагревателей, м^2
Fn=3.14*dn*ln*nn;
//Масса печи, кг
mp=(A*B*H-(H-2*bb)*(A-2*bb)*(B-2*bb)-2*a*b*bb)*rost;
 
function yd=f(tau,t)
    yd=(U^2*3.14*dn^2/(4*(2.92+0.001*t)*ln)-(tst2+131*ln(tau+1)-tos)/(bb/lst+1/ast)*Fnar-10^(-8)*5.67*En*Fn*((t+273)^4-(tos+273)^4)-10^(-8)*5.67*En*Fn*((t+273)^4-(tst2+131*ln(tau+1)+273)^4)*tau)/(nz*Cz*mz+nn*Cn*mn+Cst*mp*131*ln(tau+1));   
endfunction;
    t0=10;
    tau0=0;
    tau=0:1:150;
    y=ode(t0,tau0,tau,f);
//Вывод на экран температур:
disp('температура, с');
disp(y);
plot(tau,y);
xgrid();
xtitle('Распределение температуры во времени','Время, с','Температура, С');

function yd = f(t,T)
    yd=(U^2*3.14*dn^2/(4*(2.92+0.001*T)*ln)-(tst2+131*log(t+1)-tos)/(bb/lst+1/ast)*Fnar-10^(-8)*5.67*En*Fn*((T+273)^4-(tos+273)^4)-10^(-8)*5.67*En*Fn*((T+273)^4-(tst2+131*log(t+1)+273)^4)*t)/(nz*Cz*mz+nn*Cn*mn+Cst*mp*131*log(t+1));   
endfunction
 
    T0 = 10;
    t0 = 0;
    t  =0:1:150;
    T = ode(T0,t0,t,f);