What I said was that there was no convergence in the original problem after reformulation neither in Maple 12 nor in Maple 16.

The original code stripped to essentials:

restart;
L := .22:
lambda := 3.4:
rho := 2950:
cp := 1050:
a := lambda/(rho*cp):
h := 9:
T0 := 298:
sigma := 5.6697*10^(-8): epsilon:= .85:
PDE := diff(T(x, t), t) = lambda*(diff(T(x, t), x, x))/(rho*cp);
IBC := {-lambda*(D[1](T))(L, t) = h*(T(L, t)-T0)+sigma*epsilon*(T(L, t)^4-T0^4), T(0, t) = 1573, T(x, 0) = 298};
SOL := pdsolve(PDE, IBC, numeric, timestep = 100);
SOL:-plot(t = 80*3600, thickness = 3);
SOL:-animate(t = 0 .. 80*3600, frames = 200, thickness = 3);
p1 := SOL:-value();
p1(.22, 80*3600);
##That "worked", but the result is wrong in Maple 12, but fine in Maple 16.
####################################################
##And the reformulated version:
restart;
L := .22:
lambda := 3.4:
rho := 2950:
cp := 1050:
a := lambda/(rho*cp):
h := 9:
T0 := 298:
sigma := 5.6697*10^(-8): epsilon:= .85:
sys:= {diff(T(x, t), t) = lambda*(diff(T(x, t), x, x))/(rho*cp),T4(x,t)=T(x,t)^4};
IBC := {-lambda*(D[1](T))(L, t) = h*(T(L, t)-T0)+sigma*epsilon*(T4(L, t)-T0^4), T(0, t) = 1573, T(x, 0) = 298};
SOL := pdsolve(sys, IBC, numeric, timestep = 100);
SOL:-plot(t = 80*3600, thickness = 3);
##Here there is an error message in Maple 12 and 16:

Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0):
Newton iteration is not converging
#The only difference between Maple  12 and 16 is that Maple 12 says 't>0.' instead of 't>HFloat(0.0)'.

Nonlinear terms in the IBCs seem to get ignored (or complained about) in Maple 12.

I tried the following simplified version. First the obvious way, which compared to Maple 16 gives the wrong result:

restart;
PDE := diff(T(x, t), t) = diff(T(x, t), x, x);
IBC:={D[1](T)(1, t) = -T(1,t)^4, T(0, t) = 1, T(x, 0) = 2};
SOL := pdsolve(PDE, IBC, numeric,spacestep=.005);
SOL:-plot(t = 1, thickness = 3, colour = red);
SOL:-animate(t=0..1);
SOL:-value()(1,1);
#Now introduce T4(x,t) = T(x,t)^4:
restart;
PDE := diff(T(x, t), t) = diff(T(x, t), x, x);
IBC := {D[1](T)(1, t) = -T4(1, t), T(0, t) = 1, T(x, 0) = 2};
sys:={PDE,T4(x,t)=T(x,t)^4};
SOL := pdsolve(sys, IBC, numeric,spacestep=.005);
SOL:-plot(T,t = 1, thickness = 3, colour = red);
SOL:-animate(T,t=0..1);
SOL:-value(T)(1,1);
#Now the results agree with those from Maple 16.

#However, unfortunately when trying this last version on the original problem Maple 12 as well as Maple 16 complains:
Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0):
Newton iteration is not converging

Nonlinear terms in the IBCs seem to get ignored (or complained about) in Maple 12.

I tried the following simplified version. First the obvious way, which compared to Maple 16 gives the wrong result:

restart;
PDE := diff(T(x, t), t) = diff(T(x, t), x, x);
IBC:={D[1](T)(1, t) = -T(1,t)^4, T(0, t) = 1, T(x, 0) = 2};
SOL := pdsolve(PDE, IBC, numeric,spacestep=.005);
SOL:-plot(t = 1, thickness = 3, colour = red);
SOL:-animate(t=0..1);
SOL:-value()(1,1);
#Now introduce T4(x,t) = T(x,t)^4:
restart;
PDE := diff(T(x, t), t) = diff(T(x, t), x, x);
IBC := {D[1](T)(1, t) = -T4(1, t), T(0, t) = 1, T(x, 0) = 2};
sys:={PDE,T4(x,t)=T(x,t)^4};
SOL := pdsolve(sys, IBC, numeric,spacestep=.005);
SOL:-plot(T,t = 1, thickness = 3, colour = red);
SOL:-animate(T,t=0..1);
SOL:-value(T)(1,1);
#Now the results agree with those from Maple 16.

#However, unfortunately when trying this last version on the original problem Maple 12 as well as Maple 16 complains:
Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0):
Newton iteration is not converging

I find it desirable that the printed output distinguishes between different objects.

I find it unfortunate (but can easily live with it) that 'Pi' and 'pi' prints the same.
'e' and exp(1) are clearly distinguishable in print.

With typesetting set to 'Extended' as you describe the output from these two commands look the same, but they are not:
evalf(-a);
-a;
#Both are products, but the factors are not the same:

evalf(-a);
op(%);
eval(%%,a=5/6);
-a;
op(%);
eval(%%,a=5/6);

I find it desirable that the printed output distinguishes between different objects.

I find it unfortunate (but can easily live with it) that 'Pi' and 'pi' prints the same.
'e' and exp(1) are clearly distinguishable in print.

With typesetting set to 'Extended' as you describe the output from these two commands look the same, but they are not:
evalf(-a);
-a;
#Both are products, but the factors are not the same:

evalf(-a);
op(%);
eval(%%,a=5/6);
-a;
op(%);
eval(%%,a=5/6);

@AliKhan

I am not using 'assume' but 'assuming' as seen below.

restart;
W := m*Pi/b; K := 2*n*Pi/t;
A := simplify(2*int(sin(W*(x+1/2*b))*cos(K*x), x = -1/2*b .. 1/2*b)/t) assuming m::integer;
Aodd:=simplify(A) assuming m::odd;
simplify(A) assuming m::even;
#Supposed answer:
ANS:=4*m*b*t*cos(Pi*n*b/t)*sin((1/2)*m*Pi)^2/(Pi*(-m^2*t^2+4*n^2*b^2));
simplify(ANS) assuming m::even;
simplify(ANS) assuming m::odd;
#We notice that ANS has the wrong sign. To make a quick check:
#Numerical integration with the special values t=1,b=1,n=1,m=7:
evalf(eval(2*Int(sin(W*(x+1/2*b))*cos(K*x), x = -1/2*b .. 1/2*b)/t,{t=1,b=1,n=1,m=7}));
#Compare with the next two outputs:
evalf(eval(Aodd,{t=1,b=1,n=1,m=7}));
evalf(eval(ANS,{t=1,b=1,n=1,m=7}));
#So Maple is right: ANS has the wrong sign.
#The special case m*t = 2*n*b (notice 'Int' instead of 'int'):
Aspec := eval(2*Int(sin(W*(x+(1/2)*b))*cos(K*x), x = -1/2*b .. 1/2*b)/t, n=m*t/(2*b));
A1:=value(Aspec);
#Answers don't simplify easily
A1even:=simplify(A1) assuming m::even;
eval(A1even,m=2*p);
simplify(%) assuming p::integer;
simplify(A1) assuming m::odd;
#Proposed answer:
ANS1:=b*sin((1/2)*m*Pi)/t;
simplify(ANS1) assuming m::odd;
#At least the last two 'simplies' to the same thing! But somewhat nicer:
eval(A1,m=2*p+1);
simplify(%) assuming p::integer;
eval(ANS1,m=2*p+1);
simplify(%) assuming p::integer;

@AliKhan

I am not using 'assume' but 'assuming' as seen below.

restart;
W := m*Pi/b; K := 2*n*Pi/t;
A := simplify(2*int(sin(W*(x+1/2*b))*cos(K*x), x = -1/2*b .. 1/2*b)/t) assuming m::integer;
Aodd:=simplify(A) assuming m::odd;
simplify(A) assuming m::even;
#Supposed answer:
ANS:=4*m*b*t*cos(Pi*n*b/t)*sin((1/2)*m*Pi)^2/(Pi*(-m^2*t^2+4*n^2*b^2));
simplify(ANS) assuming m::even;
simplify(ANS) assuming m::odd;
#We notice that ANS has the wrong sign. To make a quick check:
#Numerical integration with the special values t=1,b=1,n=1,m=7:
evalf(eval(2*Int(sin(W*(x+1/2*b))*cos(K*x), x = -1/2*b .. 1/2*b)/t,{t=1,b=1,n=1,m=7}));
#Compare with the next two outputs:
evalf(eval(Aodd,{t=1,b=1,n=1,m=7}));
evalf(eval(ANS,{t=1,b=1,n=1,m=7}));
#So Maple is right: ANS has the wrong sign.
#The special case m*t = 2*n*b (notice 'Int' instead of 'int'):
Aspec := eval(2*Int(sin(W*(x+(1/2)*b))*cos(K*x), x = -1/2*b .. 1/2*b)/t, n=m*t/(2*b));
A1:=value(Aspec);
#Answers don't simplify easily
A1even:=simplify(A1) assuming m::even;
eval(A1even,m=2*p);
simplify(%) assuming p::integer;
simplify(A1) assuming m::odd;
#Proposed answer:
ANS1:=b*sin((1/2)*m*Pi)/t;
simplify(ANS1) assuming m::odd;
#At least the last two 'simplies' to the same thing! But somewhat nicer:
eval(A1,m=2*p+1);
simplify(%) assuming p::integer;
eval(ANS1,m=2*p+1);
simplify(%) assuming p::integer;

int(simplify(f),r); doesn't work for me in Maple 16, 15, and 12.

Even when specifying n success doesn't come often with or without simplify:

[seq(int(f,r),n=1..10)];
has~(%,int);

int(simplify(f),r); doesn't work for me in Maple 16, 15, and 12.

Even when specifying n success doesn't come often with or without simplify:

[seq(int(f,r),n=1..10)];
has~(%,int);

@pgr Using that  sqrt(2)*[0,0,1] is not simplified to [0,0,sqrt(2)] you could fake a third term term like this:

phi:=(x,y,z)->sqrt(2)*N+(1-sqrt(2))*N+2*([x,y,z]-N)/norm2([x,y,z]-N);

#Maple adds two list of the same length: [a,b]+[c,d] = [a+c,b+d] (automatic simplification, i.e. done before actual evaluation, so that also '[a,b]+[c,d]'; results in [a+c,b+d]). Also '2*[a,b]'; results in [2*a,2*b].
Similarly, '1.234*[a,b]'; results in [1.234*a, 1.234*b]. However, x*[a,b] is not in general simplified to [x*a,x*b], automatically or not, and that is likely the way it should be. The problem is rather with the procedure 'transform' in the plottools package, which ought to have caught the fact that phi(x,y,z) didn't evaluate to a list.
Vectors behave differently: x*; evaluates to.

Simplification of a sum of terms of the form x*[a,b] could be done like this:

#In your example:
evalindets( N+2*([x,y,z]-N)/norm2([x,y,z]-N), list, convert, Vector);
phi:=unapply(convert(%,list),x,y,z);

#In general you can define `simplify/list` like this:

`simplify/list`:=proc(u) convert(evalindets(u,list,convert,Vector),list) end proc;

#Then you would have the desired effect from

simplify( x*[a,b],list);
simplify( sqrt(2)*[a,b]+y*[a,b],list);

PS. Here is another way of 'faking' a 3rd term:

phi:=(x,y,z)->if not type([x,y,z],list(realcons)) then
            [1,2,3]
          else
            N+2*([x,y,z]-N)/norm2([x,y,z]-N)
          end if;

The effect of this is to make the count of dimensions correct. Instead of [1,2,3] you could have any list with 3 elements, or indeed any object with 3 operands, like e.g. q(k,i,j) (q unassigned), since nops(q(k,i,j)) =3.

@pgr Using that  sqrt(2)*[0,0,1] is not simplified to [0,0,sqrt(2)] you could fake a third term term like this:

phi:=(x,y,z)->sqrt(2)*N+(1-sqrt(2))*N+2*([x,y,z]-N)/norm2([x,y,z]-N);

#Maple adds two list of the same length: [a,b]+[c,d] = [a+c,b+d] (automatic simplification, i.e. done before actual evaluation, so that also '[a,b]+[c,d]'; results in [a+c,b+d]). Also '2*[a,b]'; results in [2*a,2*b].
Similarly, '1.234*[a,b]'; results in [1.234*a, 1.234*b]. However, x*[a,b] is not in general simplified to [x*a,x*b], automatically or not, and that is likely the way it should be. The problem is rather with the procedure 'transform' in the plottools package, which ought to have caught the fact that phi(x,y,z) didn't evaluate to a list.
Vectors behave differently: x*; evaluates to.

Simplification of a sum of terms of the form x*[a,b] could be done like this:

#In your example:
evalindets( N+2*([x,y,z]-N)/norm2([x,y,z]-N), list, convert, Vector);
phi:=unapply(convert(%,list),x,y,z);

#In general you can define `simplify/list` like this:

`simplify/list`:=proc(u) convert(evalindets(u,list,convert,Vector),list) end proc;

#Then you would have the desired effect from

simplify( x*[a,b],list);
simplify( sqrt(2)*[a,b]+y*[a,b],list);

PS. Here is another way of 'faking' a 3rd term:

phi:=(x,y,z)->if not type([x,y,z],list(realcons)) then
            [1,2,3]
          else
            N+2*([x,y,z]-N)/norm2([x,y,z]-N)
          end if;

The effect of this is to make the count of dimensions correct. Instead of [1,2,3] you could have any list with 3 elements, or indeed any object with 3 operands, like e.g. q(k,i,j) (q unassigned), since nops(q(k,i,j)) =3.

@Axel Vogt You may find it strange, but it works clearly the way it was intended to work. Giving the default value makes the argument optional.

Also according to specifications:

g:=proc(n::posint:=Pi) return n+1; end proc;
g(k);
g();
#If you want an error from g(k); you could close with a dollar sign:
g:=proc(n::posint:=Pi,$) return n+1; end proc;
g(k);
g();

@Axel Vogt You may find it strange, but it works clearly the way it was intended to work. Giving the default value makes the argument optional.

Also according to specifications:

g:=proc(n::posint:=Pi) return n+1; end proc;
g(k);
g();
#If you want an error from g(k); you could close with a dollar sign:
g:=proc(n::posint:=Pi,$) return n+1; end proc;
g(k);
g();

You probably meant plot([seq(Cplot( beta,lambda), beta = 0 .. 15, 5)], lambda = 0 .. 60).

But besides that it doesn't work right. It always produces the expression after 'else'.

It should be noted that the correct results can easily be made with piecewise:

lambda0 := 0.22e-1*beta^2+(7/1000000000000000)*beta+5.6;
f:=piecewise(lambda < lambda0,0.1e-2,((1/2)*lambda+(1/2)*(-1)*0.22e-1*beta^2+(1/2)*(-5.6))*exp((-1)*.17*lambda));
plot([seq(f,beta=0..30,5)],lambda=0..60);

You probably meant plot([seq(Cplot( beta,lambda), beta = 0 .. 15, 5)], lambda = 0 .. 60).

But besides that it doesn't work right. It always produces the expression after 'else'.

It should be noted that the correct results can easily be made with piecewise:

lambda0 := 0.22e-1*beta^2+(7/1000000000000000)*beta+5.6;
f:=piecewise(lambda < lambda0,0.1e-2,((1/2)*lambda+(1/2)*(-1)*0.22e-1*beta^2+(1/2)*(-5.6))*exp((-1)*.17*lambda));
plot([seq(f,beta=0..30,5)],lambda=0..60);