@Kitonum 

In case that I have more than one list the "do" loop does not seem to work:

L[1]:=[x1,x2,x3,x4]:
L[2]:=[y1,y2,y3,y4]
......

for j from 1 to 2 do
K[j]:=(i::integer) -> `if`(irem(i,4)<>0, L[j][irem(i,4)], L[j][4]);
od;

Any idea to make this procedure working for any L[i]?

@Kitonum Thank you a lot!

@Carl Love 

For fixed t=T every time I took the solutions of h(x,T) and u(x,T) and made from these polynomial functions (with cubic splines). If the result was analytic (no errors) the above curves should be zero. Above we see that as the t increases the orange curve moves away from zero (x-axis) while the blue remains the same (if not closer). As I mention I dont know if this is a reliable measure of accuracy of the whole solution but the evergoing growth of the orange curve is a sign that something goes very wrong.

@Preben Alsholm 

I dont know if it is a measure of accuracy but I did the following:

I fix t=1 and I found the values of h(x,1) and u(x,1), then with Splines of order 3 I create the functions hspline(x) and uspline(x). Then I checked if these functions are solutions of the PDE system. The orange is for Maple the blue for Mathematica.

For hspline(x):

For uspline(x):

 I dont think that anyone can come to any conclusion just by that but is a start..

EDIT

The same for t=4

hspline(x):