@janhardo You have 

uxx := diff(u[0](x), x $ 2);

but u[0] is already a function of x, so it should be

uxx := diff(u[0], x $ 2);

This is the basic way to do it, with everything in a big loop.

@salim-barzani As I said before, 

-diff(u[i], x $ 2)*I + A[i]

should have some t dependence or you can't expect to correctly take a Laplce transform wrt t. So you need to put the dependence of t in here (or just have it in from the beginning).

@salim-barzani I didn't see where you divided by 2. P[1] looks the same in both versions?

Download P1.mw

int(convert(1/A, Heaviside), x=-1..1);

@salim-barzani Your T[n] in DrD.mw has conjugate(u[k]) outside the summation, so it just stays as k.

Interesting. Nested radicals frequently give me problems. It is interesting that solve thinks they are the same.

expr:=sqrt(x+2*sqrt(x-1));
expr2:=sqrt(x-1)+1;
solve(expr2-expr,x);

gives x.

Probably you can use polygon from the plottools package for this, but if you upload your worksheet using the green up-arrow in the Mapleprimes editor, someone will help you.

@salim-barzani You didn't attempt to use Diff - all your Q's are zero because your B1 are zero.

@salim-barzani Does Diff work for you?

Download Diff.mw

@salim-barzani I didn't understand that. Can you give a very simple example of where conjugate doesn't allow you to do something which you want to do?

@salim-barzani You say "for diff(conjugate(u(x,t)),x) i use v[i,x]", but diff(conjugate(u(x,t)),x) is a function of t but v[i,x] is not. So for example the Laplace transform wrt t of v[i,x] will not be correct. If you fix these things up, I'll take another look.

I don't exactly understand what you want, but here are some remarks:

1. The calculation needs to be done iteratively (I think the paper uses recursively where I would use iteratively). So you need to work out all the n=0 things before the n=1 and so on. You can have one big loop but you cannot work out, say, A[0]..A[4], before you work out u[0](x,t)..u[4](x,t).

2. For debugging at least I would not use "declare"; it hides what is really happening. For example u[0](x,t) will appear as u[0], but so will u[0]. When you define B[0], you use u[0] but I think this has to be entered as u[0](x,t).

3. You ask about mixing u[i] with u[i,x]. It seems you should never have u[i.x] - only ever one subscript. The display of derivatives as subscripts doesn't mean they are actually subscripts (see point 2). The paper says u(0,x) is an initial condition, but that should be u(x,0) since you have functions of (x,t). So you can have u[2](x,0) but never anything like v[0,x] (which should be diff(v[0](x,t),x) or perhaps it means an initial condition v(x,0) but for which there presumably is a missing subscript.)

4. It is bad practice to use both u(x,t) and u[0](x,t) since they are not independent of each other; instead use U(x,t) and u[0](x,t). It may work if you are not going to assign to u, but things are not what they seem - try u[0](x,t):=g; then eval(u);

@Christopher2222 Well, maybe the symbols were generated by some other calculation which might or might have made them co-linear, but I agree that getting the assumptions right makes using the geometry package frustrating. The hints about what to assume are not always as clear as in this case.

@rockyicer  You can make your permutation matrix for the permutation [1,2,3,4] -> [2,1,4,3] using.

P := IdentityMatrix(4)[[2, 1, 4, 3]]

or you can just do the permutations of the rows and columns by

Y4_R := Y4[[2, 1, 4, 3], [2, 1, 4, 3]];

@acer Nice. I thought I had tried that, but I guess not.