Ошибка деление на ноль при моделировании движения ракеты относительно цели

d=0.25; %диаметр (м)
l=4.5; %длина (м)
m=360; %масса (кг)
tdv=22000; %тяга двигателя (кг)
sr=0.05; %площадь рулей (м^2) 
cy=4.8; %коэф. аэродин. силы
cx=0.41; %коэф. лоб. сопр.
la=0.5; %плечо аэродин. (м)
lr=2.25; %плечо руля (м)
h=5000; %высота полета (м)
g=9.81; %ускор. своб. падения
Kiz=10; %стат. коэф. передачи измерительного звена
Kzvk=1; %коэф. передачи звена выроботки команд
Koc=0.12; %коэф ОС
Km=0; %метод моделир. 1-погони, 0-прямое навед.
%расчет плотности воздуха по стандарту СА-81
ro0=1.225; %плотность воздуха на уровне моря
t0=288.15; %температура на уровне моря
t=t0-0.0065*h; %температура
ro=ro0*(t/t0)^4.26; %плотность воздуха
%расчет коэффициентов линеаризованной модели ЛА
s=pi*d^2/4; %площадь миделя
v=sqrt(2*tdv/cx/ro/s); %скорость ЛА
sn=ro*v^2/2; %динамическое давление
k1=s*cy*sn; %коэф. пропорц. аэродин. силы
k2=sr*sn*lr; %коэф. пропорц. актив. вр. мом.
k3=0.01*k2; %коэф. пропорц. демпфир. мом.
k4=k1*la; %коэф. пропорц. стабилиз. мом.
Iy=m*l^2/12; %мом. инерции
tv=m*v/(k1+tdv); %пост. вр. ЛА
ad=k2/Iy; %коэф. числителя
w2=(k3+k4*tv)/Iy/tv; %кв. собств. част. колебаний ЛА
d=(Iy+k3*tv)/(2*sqrt(Iy*tv*(k3+k4*tv)));%коэф. демпфир.
dlt=2*d*w2^0.5; %коэф. при первой производн.
Vtag=900 %скорость самолета;
Vp=2Vp %скорость ракеты;
N0=3 %коэффициент пропорциональности метода наведения