We can do the equivalent thing with fdiff. Here's your fdiff command:

fdiff(n*Zeta(n, 3), [n], n = 3);

We wrap any forced-integer argument with round:

fdiff(n*Zeta(round(n), 3), [n], n = 3);
      -0.3740436824

This works regardless of whether `diff/Zeta` has been modified as I suggested above because fdiff doesn't use diff anyway.

There is only 1 derivative with respect to x in your "PDE" system, and it's a first derivative. Thus 1 and only 1 condition with x fixed is allowed.

IMO, your system isn't a PDE system in the usual sense of "PDE" because each equation only has derivatives wrt a single variable. However, pdsolve(..., numeric) may solve it anyway, once you give reasonable boundary conditions. I'm curious to see. If not, then you'll need to use Preben's Heaviside workaround from the Q&A that you just linked to.

@Mariusz Iwaniuk It doesn't make sense to me to differentiate Zeta(n,z) wrt n because n is the order of the derivative. Do you mean a discrete derivative Delta(f, n) = f(n+1) - f(n)? If that's what you mean, it could be easily coded into Maple.

@Christian Wolinski The problem is that, unlike a procedure, the code of a module's "body" is not stored with the module. Indeed, it's not stored anywhere at all (other than possibly in a source-code document); it gets marked for garbage collection immediately after it's executed. All that's stored of a module is its header (as I show in my Answer below) and the values of its locals and exports at the time that the storage occurred. It's useful to keep that in mind when writing a module. Consequently, the vast majority of modules that I write either have no body at all, or the body is the single statement ModuleLoad() (whose only purpose to be in the body is to "kickstart" the module the very first time that it's run, which is useful for debugging).

@Adam Ledger 

See Donald E Knuth's The Art of Computer Programming. Specifically, half of Volume 2 is devoted to random number generators. If you were going to prison for a life sentence, and you could take just one book, this should be the one that you take.

@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 &.