Pde resolution with temperature-dependent material properties

Hi All,

I am tryng to solve a pde equation with temperature-dependent material properties and I was wondering if what I implemented was right. I have not managed to check if, for each time-step, Matlab calculates material properties using the current solution of temperature and the polynomial functions prescribed. Could you please give me a hint?

Here the main function is:

function [c,f,s] = pdex1pde(x,t,u,dudx,Tmelt,molarMass)

% Temperature-dependent material props

[ep,Kfr,cpfr,rhofr]=Kcprho(u,molarMass); %this is the function containing the polynomials for the material properties

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diffr = Kfr/(rhofr*cpfr); % frozen thermal diffusivity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c = 1/diffr; % solid phase thermal diffusivity;

f = dudx;

s = 0

end

Boundary conditions:

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t,sigma,Tenv,Boltz,molarMass)

[ep,Kfr,cpfr,rhofr] = Kcprho(ul,molarMass);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diffr = Kfr/(rhofr*cpfr); % frozen thermal diffusivity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pl = (sigma)-(Boltz*ep*(ul^4-Tenv^4));

ql = Kfr;

pr = 0;

qr = 1;

end

Initial condition:

function u0 = pdex1ic(x,Tcl)

u0 = Tcl;

end

Material property file:

function [ep,Kfr,cpfr,rhofr] = Kcprho(u,molarMass)

% Temperature-dependent material props for W

ep=0.3; %constant emissivity

%thermal conductivity in the range [300 K < u < 3695 K]

Kfr = 149.441-45.466*10^(-3)*u+13.193*10^(-6)*u^2-1.484*10^(-9)*u^3+(3.866*10^(6))/u^2; % [W/m K]

%specific heat capacity in the range [300 K < u < 3695 K]

cpfr = (21.868372+8.068661*10^(-3)*u-3.756196*10^(-6)*u^2+1.075862*10^(-9)*u^3+(1.406637*10^(-4))/u^2)/molarMass; % [J/kg K]

%mass density in the range [300 K < u < 3695 K]

rhofr = 1000*(19.25-2.66207*10^(-4)*(u-273.15)-3.0595*10^(-9)*(u-273.15)^2-9.5185*10^(-12)*(u-273.15)^3); % [kg/m^3]

end