I don't have Ralph's book and without it I cannot tell what his code attempts to do. It would have helped if you have included a statement of the mathematical problem.  In the absence of that, I have tried to guess what the problem is, and have solved it in Maple. See if that is useful to you.

But if you want to learn the details of the code in your book, I suggest that you try to apply it first to the simple problem that I stated earlier.

Download mw.mw

As Polya once wrote, for any problem that you cannot solve, there is an easier problem that you cannot solve.  Try that one first!

So I suggest that you try applying your method to solving the much simpler problem

diff(u(x,t),t) = alpha*diff(u(x,t),x,x) + 1;

with zero initial condition and zero boundary conditions.

Can you do that?

Replace your 3rd for-loop with

for i from 1 to k do
  C[i] := B[i] -~ Am[i];
end do;

The -~ combnation says that we want to subtract the scalar Am[i] from each entry of the vector B[i].

Solving x.A=b is equivalent to solving A'.x' = b', where a prime indicates the transpose.  Here is an example.

Download mw.mw

Expr := (D@@2)(T)(0) + :-O(1) + 1/2*sin(1/2*Pi*T(x))^2 + 1/12*sin(1/2*Pi*T(x))*Pi*diff(T(x), x)*cos(1/2*Pi*T(x)) + :-O(2/3*(3/2*Pi^2*diff(T(x), x)*cos(1/2*Pi*T(x))^2*diff(T(x), x, x) - Pi^3*diff(T(x), x)^3*cos(1/2*Pi*T(x))*sin(1/2*Pi*T(x)) + sin(1/2*Pi*T(x))*Pi*diff(T(x), x, x, x)*cos(1/2*Pi*T(x)) - 3/2*sin(1/2*Pi*T(x))^2*Pi^2*diff(T(x), x, x)*diff(T(x), x))/Pi^2);
convert(Expr, polynom);

Download mw.mw

Note added later:

If you prefer, you may replace the display and seq in the last few lines with plots:-animate, as in

track_points := NULL:
nframes := 30:
animate(frame, [c], c=0..2.5, frames=nframes, scaling=constrained);

doit := proc(n)
  combinat:-choose({$2..n}, 3);
  select(x->hastype(x,odd), %);
  nops(%);
end proc:

Then doit(6) returns 9.

If you would rather see the selected subsets, comment out the nops(%) line.

@MapleUser2017 

Download mw.mw

Download mw.mw

I don't understand your code for making a cylinder.  I would make a cylinder with the command in the plottools package, like this:

with(plots):
with(plottools):
display(cylinder(strips=60));

The "strips" option specifies the number of strips that make up the cylinder's cirved surface.  Change as needed. Other parameters may be specified as well.  See the help page.

You may translate and rotate the resulting cylinder through the translate and rotate commands provided in the plottools package.

The interchange of z and tau does not affect the innermost integral, therefore we need to concern ourselves with the two outer integrals.  In the picture below I have presented the calculations and the final result without words.  See if you can make sense of it.

This doesn't make much sense but it does the job:

simplify(v);

Added later: A couple of funky ways of acheiving the same result:

v + <0,0>;
v * sin(Pi/2); 

 

Your integrand blows up at x = −d/c.  The expression for the evaluated integral varies, depending on whether that singularity lies inside or outside the interval (a,b).  You can let Maple account for all possibilities by doing

int(1/(c*x + d), x = a .. b, AllSolutions); 

Or you may prefer to set the conditions yourself.  For instance, if −d/c is outside and to the left of the interval (a,b) and a < b, then we do

int(1/(c*x + d), x = a .. b) assuming a > -d/c, b > a;

But if −d/c  is inside the interval (a,b) (and a < b) then we do

int(1/(c*x + d)^2, x = a .. b) assuming real, a < -d/c, b > -d/c;

Your second integral can be handled in the same way.

There are several mathematical and programming errors in what you have shown.  Rather than pointing them one by one, I will show you a corrected worksheet.  Study it and ask if you need explanations.  The equation models the cooling of a disk of radius a (in 2D) as it loses heat through the boundary to its surrouindings.

>  restart; >  pde := diff(T(r,t), t) - 1/r*diff(r*diff(T(r,t),r), r) = Q; >  ibc := { D[1](T)(0,t) = 0,          D[1](T)(a,t) = -k*(T(a,t)-T__amb),          T(r,0) = 1 }; >  params := Q = 0, a = 1, k = 1, T__amb = 0; >  eval(pde, {params}); eval(ibc, {params}); pdsol := pdsolve(%%, %, numeric, time=t, range=0..1); The graph of the function : >  pdsol:-plot3d(t=0..2, style=patchcontour, shading="z"); Animation of the disk's tenperature versus time: >  pdsol:-value(output=listprocedure); myT := rhs(%[3]); >  plots:-animate(plot3d, [[r*cos(s), r*sin(s), myT(r,t)], r=0..1, s=-Pi..Pi],     t=0..2, style=surface, axes=framed, shading="z", frames=50);

Download heat-equation-on-a-disk.mw

Download mw.mw