4th Order Runge-Kutta NaN and NaNi Outputs

I am attempting to solve a set of coupled nonlinear differential eqns using 4th order Runge-Kutta method. r and T are unknown variables dependent on time. Values for the constants, intial conditions, and P function are in the code. When I run the code I get

Warning: Imaginary parts of complex X and/or Y arguments ignored. 
> In RungaKutta_WaterVent_20230319 (line 40) 
Warning: Imaginary parts of complex X and/or Y arguments ignored. 
> In RungaKutta_WaterVent_20230319 (line 45) 

and when I look at my values for r and T, they are NaN and NaNi. My code is below, initially I did not use dot operators, I have added them in and have seen no change. I know the equations are solvable and result in logrithmic graphs, but I cannot seem to find values using this code. I also been seeing ode45 used as a checker for these types of problems, any information on them and how I can apply it to my problem would be helpful, thank you!

%------------------ constants
p = 0.92; % [g/cm^3] ice density
C = 1.9; % [J/gK] (=0.45 cal/gK) average specific heat (at constant volume) of ice between 250 and 180 K
L = 2.4*10^3; % [J/g] average heat of sublimation of ice over this T range
Te = 280; % [K] T for the earth approximated as a 280 K blackbody
Ts = 5800; % [K] T for the sun approximated as a 5800 K blackbody
Oe = 1.4*pi; % [sr, steradian] solid angle subtended by earth at the 329 km altitude particle (taking into account the infrared opacity of the atmosphere above the hard Earth surface)
Os = 6.8*10^(-5); % [sr] solid angle subtended by the sun at the particle
s = 5.7*10^(-12); % [J/scm^2K^4] Stefan-Boltzmann constant
ke = 1.3*10^3; % [1/cm] from Mie calcs
kp = 1.3*10^3; % [1/cm]
ks = 0.13*10^3; % [1/cm]
E = (3*L)/C; % constant
y = ((3*s)/(4*pi*p*C))*((ke*Oe*(Te^4))+(ks*Os*(Ts^4))); % constant
d = (3*s*kp)/(p*C); % constant
%-----------------------------

% 4th order runga kutta method
h = 0.01;  % set the step size
t = 0:h:100;  % set the interval of t
x = zeros(2,length(t));
x(:,1) = [0.00003; 250];% set the intial value for x
n = length(x)-1;
f = @(t,x) myderiv(t,x); %insert function to be solved
for ii = 1:n
    K1 = f(t(ii),x(:,ii));
    K2 = f(t(ii)+.5.*h,x(:,ii)+K1.*.5.*h);
    K3 = f(t(ii)+.5.*h,x(:,ii)+K2.*.5.*h);
    K4 = f(t(ii)+h,x(:,ii)+K3.*h);
    x(:,ii+1) = x(:,ii)+ (h/6).*(K1+2.*K2+2.*K3+K4); 
    
end

r = x(1,:);
T = x(2,:);

figure(1)
yyaxis left
plot(t,T)
xlabel('$t [s] $',Interpreter='latex')
ylabel('$ Temperature [K] $',Interpreter='latex')
axis([0 100 170 205])
yyaxis right
plot(t,r)
ylabel('$ Radius [cm] $',Interpreter='latex')
axis([0 100 0 0.000035])
title('Runga Kutta Water Vent')



% function of ODE frm water vent
function dxdt = myderiv(t,x)
%-------------------------------constants
p = 0.92; % [g/cm^3] ice density
C = 1.9; % [J/gK] (=0.45 cal/gK) average specific heat (at constant volume) of ice between 250 and 180 K
L = 2.4*10^3; % [J/g] average heat of sublimation of ice over this T range
Te = 280; % [K] T for the earth approximated as a 280 K blackbody
Ts = 5800; % [K] T for the sun approximated as a 5800 K blackbody
Oe = 1.4*pi; % [sr, steradian] solid angle subtended by earth at the 329 km altitude particle (taking into account the infrared opacity of the atmosphere above the hard Earth surface)
Os = 6.8*10^(-5); % [sr] solid angle subtended by the sun at the particle
%Op = 4*pi; % [sr] solid angle into which the particle radiates (isotropically)
%Ti = 250; % [K] intial temperature
%ri = 3*10^(-5); % [cm] initial radius 0.3 microns
s = 5.7*10^(-12); % [J/scm^2K^4] Stefan-Boltzmann constant
ke = 1.3*10^3; % [1/cm] from Mie calcs
kp = 1.3*10^3; % [1/cm]
%kpt = 1.3*10^3; % [1/cm] @ 180 K
ks = 0.13*10^3; % [1/cm]
E = (3*L)/C; % constant
y = ((3*s)/(4*pi*p*C))*((ke*Oe*(Te^4))+(ks*Os*(Ts^4))); % constant
d = (3*s*kp)/(p*C); % constant
%-------------------------------
r = x(1);
T = x(2);
P = (2.4.*10.^10).*exp(1).^(-6110./T);
drdt = (0.27.*P)./(T.^(1.2));
dTdt = ((y.*r)-(d.*r.*T.^4)+(E.*drdt))./r;
dxdt = [drdt;dTdt];
end

Answer 1

[Torsten]

There are several possible reasons that you get NaN values.

r could become 0.

T could become negative such that sqrt(T) gives imaginary values.

You must test this by outputting r, T, drdt and dTdt in your function.

For a start, I'd suggest you use a well-tested integrator to solve your equations, namely ode45.

This will ensure that the reason for failure is due to your equations, not due to the numerical Runge-Kutta method you programmed.

[CJ]

I looked at some ode45 methods to solve for my equations, my code is below. I am still getting NaN and NaNi outputs. I believe this means the error is in my equations but I have double checked them all as well as my constants and they all appear correct. Could there be another reason for this error? I know these equations lead to a solution and their graph is this:

clear all; close all; clc;

tic;
tspan = 0:0.01:100;
r0 = 0.00003;
T0 = 250;
x0 = [r0; T0];
[t,x] = ode45(@(t,x)myderiv(t,x),tspan,x0);
rsol = x(:,1);
Tsol = x(:,2);
plot(t,Tsol)

function dxdt = myderiv(t,x)
p = 0.92; % [g/cm^3] ice density
C = 1.9; % [J/gK] (=0.45 cal/gK) average specific heat (at constant volume) of ice between 250 and 180 K
L = 2.4*10^3; % [J/g] average heat of sublimation of ice over this T range
Te = 280; % [K] T for the earth approximated as a 280 K blackbody
Ts = 5800; % [K] T for the sun approximated as a 5800 K blackbody
Oe = 1.4*pi; % [sr, steradian] solid angle subtended by earth at the 329 km altitude particle (taking into account the infrared opacity of the atmosphere above the hard Earth surface)
Os = 6.8*10^(-5); % [sr] solid angle subtended by the sun at the particle
%Op = 4*pi; % [sr] solid angle into which the particle radiates (isotropically)
%Ti = 250; % [K] intial temperature
%ri = 3*10^(-5); % [cm] initial radius 0.3 microns
s = 5.7*10^(-12); % [J/scm^2K^4] Stefan-Boltzmann constant
ke = 1.3*10^3; % [1/cm] from Mie calcs
kp = 1.3*10^3; % [1/cm]
%kpt = 1.3*10^3; % [1/cm] @ 180 K
ks = 0.13*10^3; % [1/cm]

E = (3*L)/C; % constant
y = ((3*s)/(4*pi*p*C))*((ke*Oe*(Te^4))+(ks*Os*(Ts^4))); % constant
d = (3*s*kp)/(p*C); % constant
r = x(1);
T = x(2);
P = (2.4*10^10)*exp(-6110/T);
drdt = (0.27*P)/(T.^1.2);
dTdt = ((y*r)-(d*r*T.^4)+(E*drdt))/r;
dxdt = [drdt;dTdt];
end

[Torsten]

There must be something wrong with your equations or an incompatibility of the units you use.

Using Octave gives a solution, but this solution is not reasonable.

Temperature rises immediately up to 2500 K and then goes back and settles at the initial temperature, radius explodes to about 2000 cm and also settles at this value.