ode45 is returning NaN while using time dependent equations

Main code defines some parameters and linear approximations for engine shutoff and stage seperation.

clear all
close all
clc
tmax=3*24*60*60;%max of three days (needed a max for the tspan in  ode45)
tvec=[0 tmax];
r_moon_earth=385000*1000; %meters
r_cape=6373.3*1000; %radius at latitude, meters
om_earth=360*pi/(24*60*60*180); %radians per second
v_cape=r_cape*om_earth; %initial speed at launch pad in theta direction
lg_orbit=36500*1000; %mean orbit of lunar gateway, meters
target=r_moon_earth-lg_orbit; %target height, meters (lower of two options since its in circular orbit, saves fuel)
z0=[0 0 0 om_earth];
intpos=-5*pi/180; %guess and check
t_shutoff=300000000; %guess and check
G=6.673e-11; %Nm^2/kg^2
m_earth=5.972e24; %kg
m_falcon=1420788; %kg
m_dragon=9525; %kg
m_fuel=399161; %kg
m_moon=7.34767309e22;
om_moon=1.022; %radians per second
r_moon_earth=385000*1000; %meters
f=22819*1000; %thrust of falcon heavy, newtons
Mt=m_dragon+m_fuel+m_falcon;
Ml=m_falcon+m_fuel;
tt=linspace(0,tmax,51);
F=f*(-atan(tt-t_shutoff)+1.568486298)/(3); %approximation of thurst shut off
Tr=F.*sin(90-(10e-30).^(1./(tt)).*90); %approximation of thrust direction change
Tth=F.*cos(90-(10e-30).^(1./(tt)).*90);
M=Mt-Ml*tt.^(1/4000);  %approximation for mass loss of rocket parts and fuel, kg
%equations of motion are only dependent on external forces  (thrust and
%gravity) there are no constraints on  the  motion of the rocket

[t,y]=ode45(@(t,y)My_ODE_Function(t,y,tt,Tr,M,intpos,om_moon,r_moon_earth,G,m_moon,m_earth,Tth),tvec,z0);

The ode function code creates the equations of motion using f=ma approach dmoon is to get distance of moon form rocket and find its effect

function [ydot] = My_ODE_Function(t,y,tt,Tr,M,intpos,om_moon,r_moon_earth,G,m_moon,m_earth,Tth)
a=y(3)+intpos-om_moon*tt;
dmoon=sqrt((r_moon_earth-(y(1))^2+(2*r_moon_earth*sin(a/(2*r_moon_earth)).^2))); %distance between rocket and moon
Fge=(G*m_earth*M)/(y(1)^2); %gravitational for due to earth
Fgmr=((G*m_moon*M)./(dmoon.^2))*sin(y(3)); %gravitional forces due to moon (driection  dependent)
Fgmth=-((G*m_moon*M)./(dmoon.^2))*cos(y(3));
M=interp1(tt,M,t);
Tr=interp1(tt,Tr,t);
Tth=interp1(tt,Tth,t);
Fge=interp1(tt,Fge,t);
Fgmr=interp1(tt,Fgmr,t);
Fgmth=interp1(tt,Fgmth,t);
ddr=(Tr-Fge+Fgmr)./M;
ddtheta=(Tth+Fgmth)./M;
ydot=[y(2);ddr;y(4);ddtheta];
end

ode45 is returning a vector of around 1909x4 double. the first row is actual number. The rest of the rows are NaN for all columns.

Am I missing a division by zero or is it somethign to do with the use of ode45?

Thank you very much and appreciate any help.

[joe jaracz]

figured it out, changed the z0 values to [.001 .001 .001 om_earth]

[darova]

thanks for the comment