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

@Markiyan Hirnyk Your Maple output from the latex command is different from mine. Mine does not have any internal square brackets, and mine uses underscores to indicate subscripting.

Please give a more-detailed example.