It needs a reality assumption on sigma to complete the transform:

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

Note that assuming does not work with variables in functions, so you can't use i[beam](phi), sigma being internal to i[beam]. However, assume does work with variables in functions, so you could do this also:

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

In the call to SR, you need to pass f instead of f(x). In other words, change the line

SR(0, 150, f(x), 10);

to

SR(0, 150, f, 10);

I also strongly recommend that you change the definition of f from

f(x):= (1)/(5.6 sqrt(2 Pi))*exp(-((x-169)^( 2 ))/(2(5.6)^(2)));

to

f:= x-> 1/(5.6*sqrt(2*Pi))*exp(-(x-169)^2/(2*5.6^2));

The latter is easier to read and will work regardless of the version of Maple being used.

The lower and upper bounds of a seq must be set at the time the seq is executed. You tried to execute these seq statements before the value of N had been set:

sys := [seq(Stencil[1](h,m,lambda,i[2],i[1],phi,f),m=0..2*N)];
w := [seq(phi[i],i=1..2*N+1)]:

Below that you set N:= 4. This needs to be above the seq statements.

This change doesn't solve all the problems in your worksheet, but it does get you past the GenerateMatrix step.

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});

Numeric solution:

Sol:= pdsolve({PDE}, {IBCs}, numeric):
Sol:-animate(t= 0..20, frames= 200);

Edit: See my Reply below for a much better animation.

To complete the integral, it needs to know a little about a. You can add an assuming clause to the integral:

int(int(x^2+y^2, y= -sqrt(2*a*x-x^2)..sqrt(2*a*x-x^2)), x= 0..2*a) assuming a > 0;

You can replace a > 0 with a < 0, a::real, etc.

Another option is the include the option AllSolutions:

int(int(x^2+y^2, y= -sqrt(2*a*x-x^2)..sqrt(2*a*x-x^2)), x= 0..2*a, AllSolutions);

It is also possible to get an answer for complex values of a:

a:= c+d*I:
int(int(x^2+y^2, y= -sqrt(2*a*x-x^2)..sqrt(2*a*x-x^2)), x= 0..2*a) assuming c > 0, d > 0;

The read command can be used like an include. You can read in any section of code from a plaintext file consisting of complete statements.

Try

v1:= [seq((tau^2)[lambda,i], i= 1..2)];

latex(v1);

[{\tau}^{2}_{{\lambda,1}},{\tau}^{2}_{{\lambda,2}}]

The reason for the subs problem is that the integrals on the left sides of the substitutions have minus signs in front of them. It can't find the corresponding thing with the minus signs in the substitution target. So, change Subs to

Subs:= {-b[1]= -i[1], -b[2]= -i[2]};

Student:-LinearAlgebra:-NullSpace(< < 2 | -3 | 6 > >);

There are three differences between map and LinearAlgebra:-Map (henceforth, simply Map): in-place operation, the treatmeat of immutable entries, and the options that can be passed in square brackets.

In-place operation: Map(f, A) changes A. But map(f, A)---like the vast majority of Maple commands---does not change its arguments. It makes a copy of A and changes that. So, your example that purports to show that sometimes map and Map operate the same is incorrect.

restart:
with(LinearAlgebra):
A:= <2>:
B:= Map(x-> x^2, A);

A[1];

               4

B:= map(x-> x^2, A);

A[1];

            4

The treatment of immutable entries: This is what is shown correctly by your second example. If certain entries cannot be changed, then Map does not change them; map does change them, but in the copy of course. In your example, the lower triangle and the diagonal cannot be changed because the Matrix is declared to be unit lower triangular.

Options: Both Map and map have options that can be passed in square brackets next to the command. These options are different, as described on the respective help pages.

Rather than saying that print produces no output, I'd say that the output of print is NULL and that it produces printed results as a side effect.

Also, the structures in your first example (in the original Question) are properly called lists, not arrays.

L:= [1,1,1,2,3,3,4]:
[ListTools:-SearchAll(1,L)];

[ListTools:-SearchAll(2,L)];

                              [4]

seq(myfun(V1[[1..k-1,k+1..-1]]), k= 1..numelems(V1));

If you need it in a for loop, then

for k to numelems(V1) do
     myfun(V1[[1..k-1,k+1..-1]])
end do:

There is very useful version of solve called eliminate which returns the partial solutions (i.e., certain variables eliminated from the equations) and the residual equations. Here, I re-solve those residual equations. For educational purposes, I recommend that you look at the results of eliminate below.

eqs:= {3*a+4*b = 3*x, 5*a+7*b = 7*x}:
vars:= {a,b}: newvars:= {z||(1..nops(vars))}:
Subs:= newvars =~ vars/~x:
subs(Subs, solve(eliminate(eqs union Subs, vars)[2], newvars));

A table is easier. Change the line

m:= (T, x1+x2);

to

m[T]:= [T, x1+x2];

That will make m a table. When you want to use the table, type

convert(m, list)

For example,

plot(convert(m, list));