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?

Taking the max without using abs seems not quite right. So either find max and min or find the max of the absolute value:
errfunc := abs@rhs(errfunc[3]);
Furthermore, since spacestep=0.4 it seems more relevant to use
maxerr := max(seq(errfunc(i*0.4),i=0..5));
and similarly for maxerr2.

I cannot tell if the following idea would be relevant for your situation:
Give the objective function option remember. If computation exceeds the timelimit you still have the computed values available in the memory table. You could sort the entries according to size (keeping track of the indices). You may need to trap error using try ... catch. 

@Markiyan Hirnyk  The limit has nothing to do with n. It is x. I don't think this question deserves more of my time. Sorry!

@Markiyan Hirnyk 

is the result of taking the limit, now a and b are a=(2*j-1)/n*Pi/2, b=a+Pi/n.

It is really not very deep. It uses nothing but Calculus 1 (in the US sense).

tan(n*x) tends to -infinity as x -> a from the right and tends to infinity when x tends to b from the left.
Thus knowing the limits of arctan you have the result.

@Markiyan Hirnyk In the integral the integrand is continuous  and surely limits are taken in the integral resulting in the values for a and b stated, i.e. a=(2*j-1)/n*Pi/2, b=a+Pi/n.

So yes, it is that simple.

@Markiyan Hirnyk Integration by parts:

restart;
A:=Int(u(x)*diff(v(x),x),x=a..b);
                                              
IntegrationTools:-Parts(A,u(x));
                                             

This is certainly valid on any interval [a,b] on which u and v are continuously differentiable. Continuity at a and b is only required of v in the sense that the limits limit(v(x),x=a,right) and limit(v(x),x=b,left) both exist. Similarly for u.
But arctan(2*tan(n*x)) has limit Pi/2 as x->b from the left and limit -Pi/2 as x->a from the right.