@Carl Love The Maple choice of parentheses () for grouping as well as for use in function application has its drawbacks, but it does follow ordinary mathematical practice. The 2D interface has probably just added to the problems. Mathematica made another choice.

I'm not so certain that Maple should issue an error message. It has always accepted 'weird' input like sin(x)(x) or x(x).
It plots both and integrates both numerically since real or complex numbers are accepted as constant functions. Replacing x with 2 above as in sin(2)(2) or 2(2) gives us sin(2) and 2, respectively.

Belaboring some:
u:=(exp(x)-1)(1-exp(-x));
ustar:=(exp(x)-1)*(1-exp(-x));
g:=(exp(x)-1);
plot([u,g,ustar],x=0..1,color=[red,blue,green]);
eval~([u,g,ustar],x=ln(2));
map(int,[u,g,ustar],x=0..1);
evalf(%);
int(x(x),x=0..1);
evalf(Int(x(x),x=0..1));

@acer I didn't catch that one!

@ 
The warning you got usually means that the optimization method erroneously computed the gradient to be zero.
However, to test the problem I made up my own data from the known values of alpha = 0.0035, theta = 60.
Obviously this should make for a good solution, but my point was to see if it iterated or not.
It did and it did find the solution alpha = 0.0035, theta = 60.

restart;
with(Statistics):
rho[0]:=1.33;
T0:=6;
f:=proc(T,alpha,theta) evalf(rho[0]+alpha*(T-T0)*((T-T0)/theta)^5*(Int(x^5/(exp(x)-1)(1-exp(-x)) ,x = 0 .. theta/(T-T0)))) end proc;
temp:=<seq(7..25)>;
fT:=rcurry(f,.0035,60);
resistance:=fT~(temp);
NonlinearFit(f,temp,resistance,output=[parametervalues, residuals],parameterranges=[0.001..0.01,10..100],initialvalues=[0.0027,50]);

Note:
I also tried to perturb the data slightly. NonlinearFit came up with a result here too:
X := RandomVariable(Normal(0, 1e-4));
A := Sample(X, 19):
resistance2:=resistance+convert(A,Vector[column]);
NonlinearFit(f,temp,resistance2,output=[parametervalues, residuals],parameterranges=[0.001..0.01,10..100],initialvalues=[0.0027,50]);

It should be NonlinearFit (not NonLinearFit).
You should use := for assignments.
Do not use T[0] when T is in use on its own. Use T0 (or whatever) instead.
The model has an unevaluated integral. Maple cannot do the integral, not even when the upper limit is 1.
You may want to try the operator form using something like

f:=proc(T,alpha,theta)
  evalf(rho[0]+alpha*(T-T0)*((T-T0)/theta)^5*(Int(x^5/(exp(x)-1)(1-exp(-x)) ,x = 0 .. theta/(T-T0))))
end proc;
#Test of f
f(7,.0027,50);

We don't get much further without the data.

@lcrpn 
Here is a slightly shorter version:

expand(f(x));
collect(%,sin(x));
convert(%,phaseamp,x);
res:=simplify(%);
#Incidentally the argument to cos can be simplified:
w:=op([2,1,2],res);
#Getting an idea:
identify(evalf(w));
#Showing that w+3*Pi/20 = 0:
u:=w+3*Pi/20;
expand(tan(u)); #we know that -Pi/2 < w <0 so that -Pi/2 < u < Pi/2
simplify(%);

I don't know what is going on, but noticed that if a is given a concrete value from the start (instead of assuming on it) e.g. a:=1/3 or a:=1/2 there is no difference.

There are a couple of syntax problems:

A missing od (end do) after the line
W[i]:=unapply(%,x):                 #Orthonormality process

f[j](x) should be replaced by f[j] in the line
c[io][j]:=int(f[j](x)*W[io](x),x=0..1);

I strongly suggest replacing linalg[matrix] by Matrix. In fact for the code you have I would use Vector(2^J+2) instead of Matrix(2^J+2,1). However, I have only looked at the syntax of the problem as presented in the code.

To see output from the two double loops at the end set printlevel at 2:
printlevel:=2;

There is no file attached. Use the thick green arrow in the icon panel of the editor.

@mahmood180 By applying the laplace transform to the ode you get an ode for the laplace transform. I Using dsolve on that ode you get an integral that seems rather undoable.

So I only see iteration as a way to solve your new ode. That method would also solve your previous one.
What is needed is knowledge of x(t) at some interval of length a = 1/4 then you can work your way up.
This will be an exact solution.
Similar to the previous problem I'm here assuming that x(t) = 0 for t < 0 and x(0) = x0. This implies that x(t) = x0 on t = 0..a, thus we have a start.

I posted the iterative solution in a general form as an answer and have removed the remainder of this comment.

@ateaktree The upload of the file didn't succeed. I got an error 404.

@mahmood180 Certainly. If I think I can contribute in any way I shall be happy to.
MaplePrimes is the perfect forum.

@mahmood180 
There is obviously something I am not getting. It doesn't seem to me that your T1 is a solution.

restart;
ode:=diff(x(t),t)=x(t-1)+2*t;
#If x(t)=0 for t<0 and x(0)=x0 does it not follow that x(t)=t^2+x0 for 0 <= t < 1?
eval(ode,{x(t)=t^2+x0,x(t-1)=0}); #Certainly is a solution

My guess is that people would find this a more interesting problem if they were told what you are trying to do.
You don't say much about that.
Also I wonder why you are using the linalg package at all. That one was superseded by the LinearAlgebra package with the release of Maple 6 quite a while ago.

@Carl Love I think that the question was raised in a class discussing the stability of the scheme.
In the help page (pdsolve,numeric,methods) it says about the Euler method,
"For PDEs that are even order in space, stability requires a restriction of the form
                                 xord
                          k < a h    
 for some a  depending upon the problem."
So the dramatic difference was "intended" so to speak.

What do you mean by a matrix's selfadjoint matrix?
Given the matrix H, what does that selfadjoint matrix A have to satisfy? What is its relation to H supposed to be?
You can write H (uniquely) as a sum of a hermitian matrix and an antihermitian one very easily:

H = (H+H*)/2 + (H-H*)/2

But is it just (H+H*)/2 you are after?