// example25.edp from chapter3/thermal.edp
verbosity = 0;
func u0 = 10 + 90*x/6;
func k = 1.8 * (y<0.5) + 0.2;

real ue = 25, alpha = 0.25, T = 5, dt = 0.1 ;

mesh Th=square(30,5,[6*x,y]);
fespace Vh(Th,P1);
Vh u=u0, uold;

real rad=1e-8, uek=ue+273;
Vh vold, w, v=u0-ue, b;
// v is (u-ue)
problem thermradia(v,w)
    = int2d(Th)(v*w/dt + k*(dx(v) * dx(w) + dy(v) * dy(w)))
                + int1d(Th,1,3)(b*v*w)
                - int2d(Th)(vold*w/dt) + on(2,4,v=u0-ue);

for(real t=0;t<T;t+=dt){
    vold=v;
    for(int m=0;m<5;m++){
       b= alpha + rad * (v + 2*uek) * ((v+uek)^2 + uek^2);
       thermradia;
    }
}
vold=v+ue; 
plot(vold);
interactive()