@Kitonum A more interesting result is obtained with tbas=10.

@Carl Love Ordinarily I would pick the initial guess as you have, but in this case such a choice makes the iteration arrive at the fixed point in a single step.  I picked my initial guess on purpose so that it takes a few iterations to arrive at the fixed point, and thus to illustrate how the iteration works.

As to a "best choice", I don't have a useful answer other than making the obvious point that the "best choice" is the solution itself.  In this problem it would have been u[0](t,x) = 1 + x + t + x*t, which is better than u[0](t,x) = 1 + x.

You ask "how can I iterate this PDE?" but you haven't shown the PDE.  The iterative formula in your worksheet is misformed and ambiguous, making it difficult/impossible to deduce the PDE from it.

Show the PDE.

@Carl Love Yes, you are right.

    Kbt1 and mubt1 are defined as functions but used as if they were expressions.
    Fht1 is defined as an expression but used as if it were a function.

There may be more such errors.

@Paulo Baumbach:  You can make your programming job significantly easier if you write a single statement per execution group.  Ask for explanation if you don't know what an execution group means.

Tom, you are not alone.  I am not receiving notifications either.  Something has stopped working on this website.

@Rouben Rostamian  Come to think of it, that was pretty dumb. The differential equation is linear, and is therefore solvable symbolically.  Why not just plot the symbolic solution?

Download mw.mw

@vv Certainly the u(x,t)=sin(x)*exp(t) solves the backward heat equation exactly.  However it walks a very fine line.  Nearby functions generally evolve quite differently and that becomes a major issue under numerical noise.  Change t=0..1 to t=0..3 in your example to see that. We can hardly say that the result "fully agrees with the exact solution".

@Carl Love Thanks for pointing to the OP= help page.  I had not noticed this extremely helpful addition in Maple2019.

@arshl Solving your modified problem leads to an integral of the type
which is not expressible in terms of elementary functions. This means that your problem has no symbolic solution.  If this problem arises in a real application, then you may seek a numerical solution.  But for that you need to supply the numerical values of all the parameters.

@Carl Love Your approach requires the initial condition to be twice differentiable.  We can see that by replacing the initial condition u(x,0)=x*x(x-1) by u(x,0)=g(x) and then observing the resulting expression for bounds.

Thus, for instance, this approach fails with the perfectly fine initial condition u(x,0)=Heaviside(x-1/2).

If, however, the twice differentiability condition on the initial condition is met, then your approach will produce the correct solution.

That said, let me reiterate that the original statement of the PDE is flawed.  The sign of the second derivative term is wrong.  To see the problem, in your solution change
psol:= pdsolve(pde, bounds, numeric):
to
psol:= pdsolve(pde, bounds, numeric, spacestep=s);
and watch what happens to the amplitude of the solution with s=0.01 and s=0.001.

Question: The use of 'union=' in your solution is new to me. Where is it documented?

Added later:
I see that C-style assignments such as
x := 12;
x += 1;
x *= 5;
also work.  These have been on my wish-list for a long time!

@abbastalebi Here is how:

plots:-odeplot(dsol, [eta, diff(f(eta),eta)]);      # first derivative
plots:-odeplot(dsol, [eta, diff(f(eta),eta,eta)]);  # second derivative

@Mariusz Iwaniuk That's excellent. A thumb up!  I wonder whether it's possible to recover the original integral through differentiating the answer.

@tomleslie No, I didn't mean to.  I played around with the equation a bit and, among other things, changed the sign of the exponent, but then I forgot to restore it.  Mea culpa!

That said, the function exp(t*x) grows astronomically large very quickly.  I can't think of any realistic problem leading to such a term.  I like the problem better with the minus in the exponent.

@kambiz1199 A complete answer to your question depends on the details of your worksheet that you have not shown, so people are just guessing what you may have in there.

Here is a complete worksheet where for simplicity I have taken the N matrix to be 2x4, and have taken some random expressions for N[1] and N[2].  If this works for you, then you should be able to extend it as needed.

Download mw2.mw

@arshl  We apply dsolve to find the solution v[3](t) of the differential equation.  The solution depends on many parameters such as p1, R, k11, a11, ro, cp2, etc.  There is nothing special about p1.  Maple is quite happy solving differential equations involving parameters.  Here is a simple example:

de := diff(x(t),t,t) + omega^2 * x(t) = 0;
dsolve(de, x(t));

The solution x is a function of t, so we write it as x(t).  It also depends on the parameter omega, but we don't write x(t,omega).