@JAMET Did  you try making x, y, and i local like I told you to??

Your code above still works for me using the given H2 and H5.

@tomleslie The OP's code will not work if either x or y are assigned non-name values before the presented code is run. I suspect that this was the cause of the error message.

I haven't executed it, but I see nothing wrong. What makes you say that it's not working?

There is something strange, however: Why do you compute 200 iterations and only look at the 5th?

A month ago I showed you how to set up a tridiagonal stochastic matrix in one simple command. Why don't you use that here?

@mehdi jafari There's also a command VariationalCalculus:-EulerLagrange.

@Rouben Rostamian  

Whenever I install a new version of Maple, the first command that I give is ?updates. 

Thank you for the detailed explanation of the Sobolev space. Yes, the smoothing is remarkable. So, go ahead and change the sign of the second derivative; my code will work just as well.

@radaar Let me know if you understand what I wrote or if you have any further questions. Using procedures as parameters in other procedures is a tricky issue in any language. Maple makes it much easier than most.

@Rouben Rostamian  

Regarding the smoothness of the IC: Yes, of course. I didn't intend for readers to infer that my technique was generally applicable. Instead, I was reacting to what I perceived to be Tom's implying that there was some choice involved for the new IC.

I'm glad that you noticed that my code computes D[1,1](u)(x,0) from the u(x,0) given in the original bounds.
 

Regarding the original problem being flawed: I understood your previous comments about this to mean that the system couldn't represent any practical physical situation. But are you indeed saying further that it's mathematically flawed: that it has no mathematical solution, however impractical?
 

Regarding union=: This is a new feature of Maple 2019 that's documented on the help page ?assignment. Let's let OP represent almost any binary infix operator (such as +, union, and, etc.). Then the expression

X OP= Y

is equivalent to 

X:= X OP  Y

I'm just guessing, but see if the VariationalCalculus package does what you want. See help page ?VariationalCalculus

@tomleslie Isn't D[2](u)(x,0) uniquely determined by plugging t=0 into the original PDE and applying the original initial conditions? Doing so, I get D[2](u)(x,0) = x^2*(x-1), like this:

restart:

#Original partial-differential-integro equation (PDIE) problem:
f:= (x^3+t^2*x^2-t^2*x+4*t*x-2*t+1)*exp(x*t)-x+1:
pdie:= diff(u(x,t),t) + diff(u(x,t),x,x) + u(x,t) + int(u(x,tau), tau= 0..t) = f;
bounds:= {u(0,t)=0, u(1,t)=0, u(x,0)=x*(x-1)}:

#Generate IC for D[2](u)(x,0):
bounds union= {simplify@isolate}(
   convert(value(eval(convert(eval(pdie, {t= 0}), Diff), bounds)), D), 
   D[2](u)(x,0)
);

#Convert original PDIE to a true PDE (i.e., without integrals): 
pde:= diff(pdie, t);

#Solve PDE system:
psol:= pdsolve(pde, bounds, numeric):
psol:-plot3d(u(x,t), x= 0..1, t= 0..0.167); 

@Joe Riel As it's currently coded, if your &@ associates, then it must associate to the right. Maple 19 has introduced right-associative neutral operators. They begin with && and otherwise are completely the same as those begining with single &.

@Joe Riel My choice for a name for a built-in identity function is `1`.

@Joe Riel Perhaps D[1$0] would be more clear. Personally, I'll use D[].

@Christian Wolinski Your effort to make it into a worksheet is appreciated.

@abbastalebi Word format is not very helpful. Plaintext or worksheet is best. Anyway, I think that I got all the necesssary code from Christian's worksheet transcription of your Question. 

So, what does parameter beta represent? Is it meant to be input to or output from the system?

Currently, the total differential order of the system is 3, and there's 1 parameter, beta. That's 4, yet you have 6 initial and boundary conditions. Those two things (the 4 and the 6) need to be equal. You could add two more output parameters.

What is your version of Maple?