The rules for changing the index in a sum are the same as the rules for changing a variable of integration. Likewise, the rule for expanding a sum is the same as that for integration.  There is a SumTools package, but it doesn't have simple manipulations like IntegrationTools's Expand and Change. So we change the inert sum to an inert integral, perform the manipulations, then change it back to a sum.

with(IntegrationTools):
A:= Sum(a[n]*phi(n)+b[n+1]*phi(n+2),n=0..N);

Expand(subs(Sum= Int, %));

The applyop is used below because we only want to apply the Change to the second operand.

subs([m= n, Int= Sum], applyop(Change, 2, %, m= n+2));

The change-of-variables step can also be done with PDEtools:-dchange, but that doesn't handle the Expand step.

You will get some minimal information about the steps by setting

infolevel[solve]:= 5:

before the call to solve.

In the first line, change with (plots) to with(plots). That is, get rid of the space between the words.

In the last line, change odeplots to plot, and end the line with a semicolon rather than a colon.

From your keywords, I see that maybe you think unapply has something to do with this problem. Indeed, that is true. In your current code, Maple is treating the function/procedure parameters omega and F as different from the omega and F that appear in A and B. Function/procedure parameters are always local to their procedure. The command unapply constructs the function/procedure body so that this is so.

A:= < < omega^2-5, 4 > | < 4, omega^2-5 > >:
B:= < -F, 0 >:
amp:= unapply(LinearAlgebra:-Norm(A^(-1).B, 2), omega, F);

amp(0,4);

plot(amp(omega,4), omega= 0..5, discont);

You'll need to rescale each frame so that its axes are the same as the first/largest frame. Using animate/animate3d will not help.

You need to replace each occurence of e^anything with exp(anything).

Download evalc.mw

Dena,

I see the source of your error. You have R:= dsolve*(eval.... That * should not be there.

You wrote:

But, In Fredholm_stencil   I have sum(  "int ( f(x), x=n*h..(n+1)*h), n=-N..N-1).

No, you don't have x= n*h... Three times in the worksheet you wrote nh when it should've been n*h.

To combine two integrals on the same limits, use combine.

The differences between a Matrix and a two-dimensional Array whose indices begin at 1 are miniscule. I doubt that you would ever notice an efficiency difference.

To use multicore processing, do this: First, check if your procedure myfunc is eligible for multithreading by checking ?index,threadsafe . If myfunc qualifies, do

M:= Matrix(4,5, (i,j)-> (i,j));
Threads:-Map[inplace](myfunc, M);

If myfunc does not qualify for multithreading, then do

M:= Matrix(4,5, (i,j)-> (i,j));
M:= Grid:-Map(myfunc, M);

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;