I want to solve the 1-D transient heat transfer equation

Define a 1-D geometry(a line) in y-direction

dT/dt = (k/(rho*cp))*d²T/dy²

I.C. : @t=0, temperature of all elements is 25°C

B.C. : The geometry is divided into several elements.

a) the bottom element is exposed to hot roller -kdT/dt = h(T_roller-T)

b) the upper element is exposed to atmosphere -k*dT/dt = h(T-T_ambient)

But i am getting different results in the matlab file and python file
the radiation boundary condition works perfectly for both of the files but the roller in contact boundary condition produce different results. the results differ by 4-5°C when the temperature of roller is low and the difference increases to about 10-11°C when the temperature of roller increases by 100°C

matlab code

function pde()
% Constants

Ny = 51;            % Number of spatial grid points
Nt = 5000;           % Number of time steps
velocity_x = 15/60; % line speed in m/s
% Thermal properties of different layers
layer1_thickness = 0.001;  % Thickness of layer 1 (meters)
layer1_k = 2.15;            % Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 0.001;  % Thickness of layer 2 (meters)
layer2_k = layer1_k;            % Thermal conductivity of layer 2 (W/m-K)
rho_1 = 1950;           %density of layer 1 (kg/m³)
rho_2 = 1950;           %density of layer 2 (kg/m³)
cp_1 = 2100;            %specific heat capacity of layer 1 (J/kg-K)
cp_2 = 2100;            %specific heat capacity of layer 2 (J/kg-K)
T_roller = 120;         %Temperature of hot roller (°C)
h_roller = 500;         %Convective heat transfer coefficient (W/m2-K)

% layer3_thickness = 0.002;  % Thickness of layer 3 (meters)
% layer3_k = 0.3;            % Thermal conductivity of layer 3 (W/m-K)

Ly = layer1_thickness+layer2_thickness;         % Total thickness of the system (meters)
% Ambient temperature and convection properties
T_ambient = 25;    % Ambient temperature (°C)
h_air = 12;        % Convective heat transfer coefficient of air (W/m²-K)

%Discretization in space
ny1 = ceil(layer2_thickness/Ly*Ny);
ny2 = Ny-ny1;

y1 = linspace(0,layer2_thickness,ny1);
y2 = linspace(layer2_thickness,Ly,ny2);
y = [y1,y2(2:end)];

% Initialize temperature matrix
initial_temperature1 = 25;
initial_temperature2 = 25;

radiative_flux = 150e3;
absorption = 50;
performance = 80;
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);

dia_roller = 0.5;                   %diameter of hot roller (m)
contact_angle = 180;                %contact angle/wrap angle for hot_roller in °

y_1 = pi*dia_roller*contact_angle/360;              %in contact with hot roller
y_2 = 1.3;                                          %after hot roller
%y_3 = 0.3;
%y_4 = 1.0;

%p0 = y_1;
%p1 = p0+y_2;
%p2 = p1+y_3;
%p3 = p2+y_4;

t1 = y_1/velocity_x;
t2 = y_2/velocity_x;
%t3 = y_3/velocity_x;
%t4 = y_4/velocity_x;

T = t1+t2;

nt1 = ceil(t1/(t1+t2)*Nt);
nt2 = Nt-nt1;
%nt3 = ceil(t3/(t1+t2+t3+t4)*Nt);
%nt4 = Nt - nt1-nt2-nt3;

t11 = linspace(0,t1,nt1);
t12 = linspace(t1,t1+t2,nt2);
%t13 = linspace(t1+t2,t1+t2+t3,nt3);
%t14 = linspace(t1+t2+t3,t1+t2+t3+t4,nt4);


%pdepe settings
m = 0;  %for 1-D cartesian coordinates with no symmetry
phase = 1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t11);
u1 = sol(:,:,1);
phase = 2;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t12);
u2 = sol(:,:,1);
%phase = 3;
%sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t13);
%u3 = sol(:,:,1);
%phase = 4;
%sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t14);
%u4 = sol(:,:,1);

plot([t11,t12],[[u1(:,1);u2(:,1)],[u1(:,25);u2(:,25)],[u1(:,50);u2(:,50)]])
grid on
function [c f s] = pdefun(y,t,u,dudy)
  if y <= layer2_thickness
    c = rho_2*cp_2;
    f = layer2_k*dudy;
    s = 0;
  else
    c = rho_1*cp_1;
    f = layer1_k*dudy;
    s = 0;
  end
end
function u = icfun(yq)
  if phase==1
    if yq <= layer2_thickness
      u = initial_temperature1;
    else
      u = initial_temperature2;
    end

  else
    u = interp1(y,u1(end,:),yq);
  end
end
function [pl ql pr qr] = bcfun(yl,ul,yr,ur,t)
  if phase==1
    %pl = -h_roller*(ur-T_roller); %% in contact with roller at bottom element
    pl = Net_radiative_intensity; %% radiation at bottom element
    %pl = -h_air*(ul-T_ambient);
    ql = 1;   
    pr = h_air*(ul-T_ambient);
    %pr = -Net_radiative_intensity;   
    qr = 1;

  %elseif phase==2
    %pl = -h_air*(ul-T_ambient);
    %ql = 1;
    %pr = h_air*(ur-T_ambient);
    %qr = 1;

  %elseif phase==3
    %pl = -h_air*(ul-T_ambient);
    %ql = 1;
    %pr = -Net_radiative_intensity;
    %qr = 1;

  else
    pl = -h_air*(ul-T_ambient);
    ql = 1;
    pr = h_air*(ur-T_ambient);
    %pr = -Net_radiative_intensity;
    qr = 1;
  end
end
end

matlab_radiation_bottom_element
matlab_roller_contact_bottom_element

python code

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
Ny = 1000           # Number of spatial grid points
Nt = 5000         # Number of time steps
velocity_x = 15 / 60  # Line speed in m/s

# Thermal properties of different layers
layer1_thickness = 0.001  # Thickness of layer 1 (meters)
layer1_k = 2.15          # Thermal conductivity of layer 1 (W/m-K)
rho_1 = 1950              # Density of layer 1 (kg/m³)
cp_1 = 2100              # Specific heat capacity of layer 1 (J/kg-K)

#layer2_thickness = 0.002  # Thickness of layer 2 (meters)
#layer2_k = 0.05          # Thermal conductivity of layer 2 (W/m-K)
#rho_2 =150              # Density of layer 2 (kg/m³)
#cp_2 = 2100              # Specific heat capacity of layer 2 (J/kg-K)

layer2_thickness = 0.001
layer2_k = 2.15
rho_2 = rho_1
cp_2 = cp_1

T_roller = 120           # Temperature of hot roller (°C)
h_roller = 500           # Convective heat transfer coefficient (W/m²-K)

T_ambient = 25           # Ambient temperature (°C)
h_air = 12               # Convective heat transfer coefficient of air (W/m²-K)

Ly = layer1_thickness + layer2_thickness

# Discretization in space
ny1 = int(np.ceil(layer2_thickness / Ly * (Ny)))  # Number of points in layer 2
#print(ny1)
ny2 = Ny - ny1
#print(ny2)                                        # Number of points in layer 1
y1 = np.linspace(0, layer2_thickness, ny1)
y2 = np.linspace(layer2_thickness, Ly, ny2 + 1)[1:]   # Exclude first point to avoid overlap
y = np.concatenate((y1, y2))
#print(y1)
#print(y2)
#print(y)
#print(len(y))
#print(y[Ny-1])
# Initial condition function
def icfun(y, ny1, ny2):
    u_initial = np.zeros_like(y)
    u_initial[:ny1] = 25  # Temperature for layer 2
    u_initial[ny1:] = 25  # Temperature for layer 1
    #print(u_initial)
    #print(len(u_initial))
    return u_initial


# Initial condition
initial_condition = icfun(y, ny1, ny2)
#print(initial_condition)
#print(len(initial_condition))

def pdefun(t, u):
    dudt = np.zeros_like(u)
    
    # Internal nodes for layer 2
    for i in range(1, Ny-1):
        #print(y[i])
        #print(i)
        if y[i] <= layer2_thickness:  # Layer 2
            dudt[i] = layer2_k / (rho_2 * cp_2) * (u[i+1] - 2*u[i] + u[i-1]) / (y[i-1] - y[i])**2
            #print(i)
        else:  # Layer 1
            #print(i)
            dudt[i] = layer1_k / (rho_1 * cp_1) * (u[i+1] - 2*u[i] + u[i-1]) / (y[i-1] - y[i])**2
    
    # Boundary conditions
    if t <= t1:
        #dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 - h_air / (rho_2 * cp_2) * (u[0] - T_ambient) / (y[1] - y[0])
        #dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 - h_roller / (rho_2 * cp_2) * (u[0] - T_roller) / (y[1] - y[0])     ##in contact with roller at bottom element
        dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 + Net_radiative_intensity / (rho_2 * cp_2 * (y[1] - y[0]))           ##radiation at bottom element
        dudt[-1] = layer1_k / (rho_1 * cp_1) * (u[-2] - u[-1]) / (y[-1] - y[-2])**2 - h_air / (rho_1 * cp_1) * (u[-1] - T_ambient) / (y[-1] - y[-2])
        #dudt[-1] = layer1_k / (rho_1 * cp_1) * (u[-2] - u[-1]) / (y[-1] - y[-2])**2 + Net_radiative_intensity / (rho_1 * cp_1 * (y[-1] - y[-2]))
        
    else:
        dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 - h_air / (rho_2 * cp_2) * (u[0] - T_ambient) / (y[1] - y[0])
        dudt[-1] = layer1_k / (rho_1 * cp_1) * (u[-2] - u[-1]) / (y[-1] - y[-2])**2 - h_air / (rho_1 * cp_1) * (u[-1] - T_ambient) / (y[-1] - y[-2])
        #dudt[-1] = layer1_k / (rho_1 * cp_1) * (u[-2] - u[-1]) / (y[-1] - y[-2])**2 + Net_radiative_intensity / (rho_1 * cp_1 * (y[-1] - y[-2]))
        
        
    return dudt



# Solve the PDE
dia_roller = 0.5
contact_angle = 180

radiative_flux = 150e3
absorption = 50
performance = 80
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100)

y_1 = np.pi * dia_roller * contact_angle / 360
y_2 = 1.3

# Discretization in time
t1 = y_1 / velocity_x
t2 = y_2 / velocity_x
T = t1 + t2
t_eval = np.linspace(0, T, Nt)
sol = solve_ivp(pdefun, (0, T), initial_condition, method='BDF', atol=1e-4, rtol=1e-4, t_eval=t_eval)

# Extract the solution
u = sol.y.T

# Plot results
plt.plot(t_eval, u[:, 0], label='Temperature at bottom surface')
plt.plot(t_eval, u[:, ny1], label='Temperature at interface')
plt.plot(t_eval, u[:, -1], label='Temperature at upper surface')
plt.xlabel('Time (s)')
plt.ylabel('Temperature (°C)')
plt.legend()
plt.show()

python_radiation_bottom_element
python_roller_contact_bottom_element

I tried to check whether the time points are defined exactly or not and checked with radiation boundary condition which works alright. But the strange behaviour is that the curve shape or behaviour is same.

Could anyone please help or suggest how can I resolve this discrepancy?
Thanks in advance!