So my first advice would be to choose some parameters that actually are the ones in the figures. Then, as @mmcdara suggested, make some approximate solutions, as I did here. This worked, but revealed a problem with one of the boundary conditions. Perhaps with this fixed, you will not need an approx solution.

Download dsolve.mw

I'm not sure why it is so slow. Note you should always use add instead of sum (convert to 1-D notation to see that you have sum).

Here's a simple way that just uses coeff. You could modify it to add conjugate etc, but it might be more efficient to substitute these in later.

Download Adom.mw

Entering control key (command on Mac) and = in the region with the integral generates the following output, so gives part of what you want. This is also available on the right-click context menu as "evaluate and display inline".

Download 5.3-Definite-Integral.mw

I ran this system with the latest version of my code, with roughly similar times to yours (and same results).

There are three calculations in the main loop - the laplace calculations, calculating and substituting into the Adomian polynomials, and evaluating the linear part. Code profiling shows that for the non-fractional order case, most of the time is spent evaluating the linear part, much to my surprise. (I'd expected the laplace calcs would be limiting.) At the moment, I can't see a way to speed this up - you need to take derivatives and simplify - if you don't simplify at each iteration then, as @acer pointed out, the u[i] components quickly become too large.

For the fractional-order case the time is divided between the linear part and the Adomian polynomials. The Adomian procedure can perhaps be speeded up; I'll look a bit more closely at that.

Download LADtest.mw

Download polylog.mw

Perhaps there is some error in pdetest, or you have a better method of verifying the solution, but my guess is the exact solution is incorrect.

Download num-Ex1_soln.mw

To plot the exact and approx solutions together in 3d

plot3d([abs(exact_soln), abs(approx_soln)], x = -2 .. 2, t = 0 .. 1, color = [red, blue]);

or Re or Im instead of abs.

For 2d, set one of x or t to a given value and do similarly with plot.

Inverse Laplace transforms often require assumptions to work. Here assuming alpha is positive works, at least for Ex 1.

Ex1.mw

file menu, bottom right corner - labeled options.

@salim-barzani I didn't use your ugly code because you did not have the correct x,t dependence in it at the time, though you since fixed that. But I also wanted general code - there is no reason to have special P,Q,R for every case when you can go directly to A[n].

Thanks for pointing out the missing factorial. You may be suggesting something else is missing but if so I didn't understand what you meant. So I think this code is correct now. The conclusion is now the same as @acer's in the thread here, i.e., that the paper has an error.

Download Laplace_Adomiian_Decomposition_method_4.mw

The error message indicates a problem with floats. Ideally, you should do the whole thing symbolically and then substitute in values, but that takes too long, probably because the expressions get too complicated. So you can convert the float parameters to rational; run through the loop and then evalf the results. Then you will need to change the parameters and run it again. It's possible that with numerical parameters you don't need the assumptions and/or the simplify (one or both), you you could try that for faster operation.

error-problem2.mw

It means there are some assumptions involving these variables. You can avoid displaying them with Options->Display->Assumed variables and choose No Annotation (or you might prefer Phrase).

There is no BlockDiagonalMatrix command. When the output of a command has the same command in it, that is a hint that the command does not exist. However, DiagonalMatrix does what you want.

Y8 := DiagonalMatrix([Y4, Y4_R]);

I don't understand the logic of what you are trying to do well enough to answer your other questions.

pdetest_explained.mw

pdetest_fixed.mw

Download parts.mw

Edit: version that accepts the set of parts as a third argument

f := proc(p::posint, n::{posint, range(posint), list(posint)}, parts::set(nonnegint))
  local sys, letters, s, w, addition, ser, z, u, letts, _s, i, cu;
  letts := seq(cat(_s, i), i in parts);
  sys := {letters = 'Union'(letts), w = 'Set'(letters), letts =~ 'Atom'};
  addition := (map2(apply, s, {letts}) =~ parts) union {s(w) = 'Set'(s(letters))};
  Order := p + 1;
  ser := combstruct:-agfseries(sys, addition, 'unlabeled', z, [[u, s]])[w(z, u)];
  ser := collect(convert(ser, 'polynom'), [z,u], 'recursive');
  cu := coeff(ser, z, p);
  if n::posint then
    coeff(cu, u, n)
  else
    [seq(coeff(cu, u, i), i = n)]
  end if
end proc:

e.g., f(15, 100, {0, 1, 4, 6, 8, 9}); gives 129.

f(15, 90 .. 100, {0, 1, 4, 6, 8, 9}); gives [201, 187, 186, 175, 172, 159, 158, 147, 143, 133, 129]

f(15, [90, 92, 100], {0, 1, 4, 6, 8, 9}); gives [201, 186, 129]

The following is just a programatic version of how you made U. When you say "generalize to higher dimensions" I'm assuming you mean larger matrices, not things in more than two dimensions that wouldn't be matrices.

MMflatten := proc(Q::Matrix(Matrix))
  local m, n, i, j;
  m, n := op(1, Q);
  <seq(`<|>`(seq(Q[i, j], j = 1 .. n)), i = 1 .. m)>
end proc: