Just one more:

evalindets(current_marks,list(integer),x->`+`(op(x))/nops(x));

@ilods It seems to be a simple zero/zero problem occurring when t (or s) is zero, so your "dirty" solution should be OK. However, the problem can be removed by inserting a line after the definition of k. Here is the line defining k followed by the new one:

k:=eval(B,{x=c[1],y=c[2]});
k:=map(normal~,k);

@ilods It seems to be a simple zero/zero problem occurring when t (or s) is zero, so your "dirty" solution should be OK. However, the problem can be removed by inserting a line after the definition of k. Here is the line defining k followed by the new one:

k:=eval(B,{x=c[1],y=c[2]});
k:=map(normal~,k);

@LijiH Yes, and it does so right away as you have also observed. But apparently you are not satisfied with that? Why not?

@LijiH Yes, and it does so right away as you have also observed. But apparently you are not satisfied with that? Why not?

@Chrissieh18 You could upload a worksheet.

Actually I strongly suspect that you are missing a multiplication sign between r and t in the exponential function. That is also the reason for the way I formulated my answer.

@Chrissieh18 You could upload a worksheet.

Actually I strongly suspect that you are missing a multiplication sign between r and t in the exponential function. That is also the reason for the way I formulated my answer.

@ilods I don't know how you modified my code to make it eat 4 giga.

I tried myself and have uploaded the result. The time used for different a's on my machine is between 2 and 3 minutes.

MaplePrimes12-09-21.mw

##############

Sept.24: Here is yet another version, which computes a plot at the average of 1 minute.

MaplePrimes12-09-21C.mw

@ilods I don't know how you modified my code to make it eat 4 giga.

I tried myself and have uploaded the result. The time used for different a's on my machine is between 2 and 3 minutes.

MaplePrimes12-09-21.mw

##############

Sept.24: Here is yet another version, which computes a plot at the average of 1 minute.

MaplePrimes12-09-21C.mw

@ilods odeplot uses adaptive plotting when you use the refine option. However, if you don't it seems to have equal spacing between the t-values. The following example seems to support this claim.

#First without range and refine:
restart;
res:=dsolve({diff(y(t),t)=x(t)^2,diff(x(t),t)=-y(t),x(0)=0,y(0)=1},numeric):
p:=plots:-odeplot(res,[x(t),y(t)],0..3.2):
plottools:-getdata(p):
M:=%[-1];
plots:-odeplot(res,[t,t],0..3.2): #To see the t-values
plottools:-getdata(%);
Mt:=%[-1];
op(1,Mt)[1];
#The spacings are the same:
{seq(Mt[i,1]-Mt[i-1,1],i=2..op(1,Mt)[1])};
3.2/200;
#Now with range and refine:
restart;
res:=dsolve({diff(y(t),t)=x(t)^2,diff(x(t),t)=-y(t),x(0)=0,y(0)=1},numeric,range=0..3.2):
p:=plots:-odeplot(res,[x(t),y(t)],0..3.2,refine=1):
plottools:-getdata(p):
M:=%[-1];
plots:-odeplot(res,[t,t],0..3.2,refine=1):
plottools:-getdata(%);
Mt:=%[-1];
op(1,Mt)[1];
#The spacings are quite different:
evalf[4]({seq(Mt[i,1]-Mt[i-1,1],i=2..op(1,Mt)[1])});
####################################
#Finally, you also have the possibility of using output=Array as in the following:
restart;
res:=dsolve({diff(y(t),t)=x(t)^2,diff(x(t),t)=-y(t),x(0)=0,y(0)=1},
numeric,output=Array([seq(3.2/200*i,i=0..200)]));
plots:-odeplot(res,[x(t),y(t)],0..3.2);
res[1,1];
TXY:=res[2,1];
TXY[1..10,1];

@ilods odeplot uses adaptive plotting when you use the refine option. However, if you don't it seems to have equal spacing between the t-values. The following example seems to support this claim.

#First without range and refine:
restart;
res:=dsolve({diff(y(t),t)=x(t)^2,diff(x(t),t)=-y(t),x(0)=0,y(0)=1},numeric):
p:=plots:-odeplot(res,[x(t),y(t)],0..3.2):
plottools:-getdata(p):
M:=%[-1];
plots:-odeplot(res,[t,t],0..3.2): #To see the t-values
plottools:-getdata(%);
Mt:=%[-1];
op(1,Mt)[1];
#The spacings are the same:
{seq(Mt[i,1]-Mt[i-1,1],i=2..op(1,Mt)[1])};
3.2/200;
#Now with range and refine:
restart;
res:=dsolve({diff(y(t),t)=x(t)^2,diff(x(t),t)=-y(t),x(0)=0,y(0)=1},numeric,range=0..3.2):
p:=plots:-odeplot(res,[x(t),y(t)],0..3.2,refine=1):
plottools:-getdata(p):
M:=%[-1];
plots:-odeplot(res,[t,t],0..3.2,refine=1):
plottools:-getdata(%);
Mt:=%[-1];
op(1,Mt)[1];
#The spacings are quite different:
evalf[4]({seq(Mt[i,1]-Mt[i-1,1],i=2..op(1,Mt)[1])});
####################################
#Finally, you also have the possibility of using output=Array as in the following:
restart;
res:=dsolve({diff(y(t),t)=x(t)^2,diff(x(t),t)=-y(t),x(0)=0,y(0)=1},
numeric,output=Array([seq(3.2/200*i,i=0..200)]));
plots:-odeplot(res,[x(t),y(t)],0..3.2);
res[1,1];
TXY:=res[2,1];
TXY[1..10,1];

What is Eq75 and what is Z2? Is Eq75 really the same as Sol?

First: Data about Ti is missing.

But apart from that you clearly get results (with any reasonable Ti) and with h=1. So what is the requirement that should enable you to determine h, i.e. the temperature dependence of the heat transfer coefficient?
Do you mean to insert h as a function of the surface temperature, like h=f(Theta(rb,t)) where f is known?

First: Data about Ti is missing.

But apart from that you clearly get results (with any reasonable Ti) and with h=1. So what is the requirement that should enable you to determine h, i.e. the temperature dependence of the heat transfer coefficient?
Do you mean to insert h as a function of the surface temperature, like h=f(Theta(rb,t)) where f is known?

Could we have the data? So we could see what is going on?

Variables := ExcelTools:-Import("C:\\Users\\Documents\\Maple\\MAPLE - MARCH 2011\\InputVariables.xls", "InputVariables", "D5:D25");