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;

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