@casperyc asked:

First of all, I have to say that 'myfun' does call a lot of other functions and procedures and returns a single value. Does that mean I have to check every single one of them?

Yes, you have to check them all. But the list is alphabetical, so it's fairly easy.

If just to be safe, say I assume it does not qualify for multithreading and use the Grid option, is there going to be less efficient? Broadly speaking, how much difference?

The difference is huge. The Grid option starts a separate mserver on every processor and copies the entire memory of your base process to each one. Threads, on the other hand, does everything in the same memory space.

@Carl Love Here is a much smoother animation which shows that the analytic and numeric solutions are in perfect agreement and which eliminates the spurious "epi-waves" caused by the default coarse grid spacing.

restart:
PDE:= diff(u(x,t), t$2) = diff(u(x,t), x$2):
IBCs:= u(x,0)=sin(Pi/2*x)*exp(-x), D[2](u)(x,0)=0, u(0,t)=0, u(4,t)=0:
Sol:= pdsolve({PDE, IBCs}):
Sol:= subs(op([2,2,1], Sol)= n, Sol):  #Change _Z1~ to n, just for aesthetics
rhs(Sol);

Sol:= subs(infinity= 20, rhs(Sol)):
A1:= plots:-animate(plot, [Sol, x= 0..4], t= 0..8, frames= 200, thickness= 3):

Numeric solution:
Sol:= pdsolve({PDE}, {IBCs}, numeric, spacestep= .02):
A2:= Sol:-animate(t= 0..8, frames= 200, thickness= 3):
plots:-display(A1,A2);
By displaying the two animations together, we can see that they are in perfect agreement because there is only ever one line.

@Joe Riel Note that the eigenvalues of your Cx are [0, 1, 1] regardless of your choice c12 or c13.

@casperyc Okay, I understand your question now. No, there is no predefined way to do that. It seems very complicated, akin to the identify command.

@casperyc Sorry, I don't understand that any better than the original. Can you give an actual example rather than an explanation?

@Mac Dude Preben's solution does not strictly require eval but it does require evaluation before passing to inttrans[fourier]. This is because assuming does not work with variables in functions. So, using evaluation, this works:

i[beam]:= phi-> exp(-(phi-phi0)^2/2/sigma^2):
ex:= i[beam](phi):
i[beamspec]:= inttrans[fourier](ex, phi, omega) assuming sigma>0;

Please see my Answer below, which you seem to have overlooked.

@Kitonum The animation quality can be improved dramatically by including spacestep= .02 in the call to pdsolve. Note that setting higher values of numpoints does not change the discretization error. Using spacestep=.02 gives 200 points on the interval x= 0..4.

@Markiyan Hirnyk Axel's code works for me. I get the same numeric answer in 10 seconds total (mostly for the integration). I am one of the upvoters.

@Carl Love Using Preben's code, I now understand why f(1) is -x (actually x-> -x) and f(x) is -1. What was confusing me was that I thought that P(x) was always 1+x+x^2.

@Alejandro Jakubi While we cannot construct a patterned matrix of unspecified (abstract) size in Maple, all the OP wants is a procedure that constructs it for specified sizes. This is quite doable in Maple.

@xcyborg You can change the options on this plot, of course. I just made some stuff up.

plot(
     [zip(`[]`, W, M), zip(`[]`, W, B)],
     style= point, symbol= [diamond, box],
     symbolsize= 15, color= [brown, green]
);

@mahmood180 

Replace the G2 loops by the single statement:

G2:= Matrix(r, (i,j)-> `if`(i+j <= r, z[i+j+r+1], 0));

When you use a[k]^`*`, do you mean the complex conjugate of a[k]? Or is it just a notation to indicate something different from a[k]?

You say that the expected result of f(1) is -x and that the expected result of f(x) is -1. Neither of those make sense to me. Could you explain why you expect those results? Then maybe I can help you with the rest. Writing a function within a function is easy.

I just need a mathematical explanation; no knowledge of or use of Maple is needed.

@ecterrab Can you provide an explicit example showing the difference between PDEtools:-Solve and solve to which you are referring? I get identical RootOf results from both PDEtools:-Solve(Theta=F(s), s) and solve(Theta=F(s), s).

@mahmood180

You need to correct the Error, Vector index out of range. Vector z has 2*r+1 elements, but the code z[i+j+r+1] is trying to access z[2*r+2] (because i = r-1 and j = 2).