Why am i getting NaN values in my code?

I'm writing a code that predicts the surface temperature of an annulus with specific thickness running through it a hot air and outisde it is under a water bath of a specific water temperature, the code should give me the way the temperature of the surface behaves across it's length.

% Input parameters
h_air = 500; % Heat transfer coefficient of hot air [W/m^2*K]
T_air_in = 600; % Inlet temperature of hot air [K]
T_water = 376; % Temperature of surrounding water [K]
k_cylinder = 50; % Thermal conductivity of cylinder material [W/m*K]
r = 0.05; % Cylinder radius [m]
L = 0.1; % Cylinder length [m]
dx = 0.005; % Grid size [m]
dt = 0.1; % Time step [s]
nt = 3600; % Simulation time [s]
d = 0.01; % Cylinder thickness [m]
h_water = 10000; % Heat transfer coefficient of surrounding water [W/m^2*K]

% Calculate grid dimensions
nx = round(L/dx) + 1;

% Set initial and boundary conditions
T = T_water * ones(nx, 1);
T(1) = T_air_in;
T(end) = T_water;

% Calculate constants for Euler's method
alpha = k_cylinder/(dx^2);
beta = h_air*dx/k_cylinder;

% Perform time integration using Euler's method
for i = 1:nt
    T_old = T;
    for j = 2:nx-1
        % Calculate the thermal resistance of the pipe wall
        r_pipe = log((r+d)/r)/(2*pi*k_cylinder*L);
        
        % Calculate the thermal resistance of the fluid film
        r_film = 1/(h_air*2*pi*r*dx);
        
        % Calculate the heat transfer coefficient
        h = 1/(r_pipe+r_film);
        
        % Calculate the temperature gradient
        dTdr = (T_old(j+1)-T_old(j-1))/(2*dx);
        
        % Calculate the temperature at the current time step
        T(j) = T_old(j) + alpha*dt*(T_old(j+1)-2*T_old(j)+T_old(j-1)) + h*2*pi*r*dt*(T_air_in-T_old(j)) + 2*pi*r*h_water*dt*(T_water-T_old(j));
    end
end

% Plot temperature profile
x = linspace(0, L, nx);
plot(x, T);
xlabel('Position [m]');
ylabel('Temperature [K]');
title('Temperature profile of cylinder surface');

Upon running the code, you'll see that it won't show a curve nor values across the cylinder.

[Torsten]

Your update

T(j) = T_old(j) + alpha*dt*(T_old(j+1)-2*T_old(j)+T_old(j-1)) + h*2*pi*r*dt*(T_air_in-T_old(j)) + 2*pi*r*h_water*dt*(T_water-T_old(j));

cannot be correct because the unit of e.g.

alpha*dt*(T_old(j+1)-2*T_old(j)+T_old(j-1))

is

W/(m^3*K) * s * K = J/m^3

but it should be

K

Similar problems with the other terms.