I'm seeing strange behavior when I try to numerically solve the heat equation with radiative boundary conditions (heat flux proportional to 4th power of surface temperature). I set up a simple example of a plate that sees 300K on one side and 310K on the other. Two obvious flaws in the solution of the radiative example are that the slope of the curve always has the wrong sign at the left boundary (sometimes wrong on both boundaries), and that the solution runs away to a temperature higher than 310K, rather than reaching a steady state between the two end temperatures. For comparison, I also ran a similar example of convective cooling (heat flux proportional to 1st power of surface temperature), which gives reasonable results. Sorry, I couldn't figure out how to post the plots.

# radiative heating at faces of plate

# density in kg/m^3

rho := 3000:

# heat capacity in J/kg*K

Cp := 1000:

# thermal conductivity in W/m*K

K := 20:

# ambient temperatures in K

T1 := 300: T2 := 310:

# Stefan-Boltzmann in W/m^2*K^4

sigma := 5.67e-8:

# thickness in m

d := 1:

DZ := d/100: DT := 1:

PDE := rho*Cp*diff(Temp(z,t),t)=K*diff(Temp(z,t),z,z);

IBC := {Temp(z,0)=T1, D[1](Temp)(0,t)=sigma*((Temp(0,t))^4-T1^4)/K, D[1](Temp)(d,t)=-sigma*((Temp(d,t))^4-T2^4)/K};

pds := pdsolve(PDE,IBC,numeric,spacestep=DZ,timestep=DT):

p1 := pds:-plot(t=100*DT,z=0..d,numpoints=50, color=red):

p2 := pds:-plot(t=1000*DT,z=0..d,numpoints=50, color=green):

p3 := pds:-plot(t=10000*DT,z=0..d,numpoints=50, color=blue):

p4 := pds:-plot(t=100000*DT,z=0..d,numpoints=50, color=black):

plots[display]({p1,p2,p3,p4}, title=`Temperature at t=100s,1000s,10000s,100000s`);

# convective heating at faces of plate

# density in kg/m^3

rho := 3000:

# heat capacity in J/kg*K

Cp := 1000:

# thermal conductivity in W/m*K

K := 20:

# ambient temperatures in K

T1 := 300: T2 := 310:

# thickness in m

d := 1:

# heat transfer coefficient in W/m^2*K

h := 100:

DZ := d/100: DT := 1:

PDE := rho*Cp*diff(Temp(z,t),t)=K*diff(Temp(z,t),z,z);

IBC := {Temp(z,0)=T1, D[1](Temp)(0,t)=h*(Temp(0,t)-T1)/K, D[1](Temp)(d,t)=-h*(Temp(d,t)-T2)/K};

pds := pdsolve(PDE,IBC,numeric,spacestep=DZ,timestep=DT):

p1 := pds:-plot(t=100*DT,z=0..d,numpoints=50, color=red):

p2 := pds:-plot(t=1000*DT,z=0..d,numpoints=50, color=green):

p3 := pds:-plot(t=10000*DT,z=0..d,numpoints=50, color=blue):

p4 := pds:-plot(t=100000*DT,z=0..d,numpoints=50, color=black):

plots[display]({p1,p2,p3,p4},title=`Temperature at t=100s,1000s,10000s,100000s`);