According to the help page ?with:

The with command is effective only at the top level, and intended primarily for interactive use. Because with operates by using lexical scoping, it does not work inside the bodies of procedures, module definitions, or within statements. See the examples at the end of this help topic.

Bottom line: you should never use "with" inside a procedure body.  Instead, you can use "uses".  Thus:

SubsColumn: = proc (A, V, k) 
local B, m; 
uses LinearAlgebra; 
m: = ColumnDimension (A); 
B: = <A[..,1..m]|V>; 
B: = ColumnOperation (B, [k, m +1]); 
B: = DeleteColumn (B, m +1); 
B; 
end proc;

If you want the heat flux to be continuous across the boundary, you need not only the thermal diffusivity r but the thermal conductivity, which is the ratio (heat flux)/(temperature gradient).  By just using the heat equation as you did, you are in effect assuming the thermal conductivity is equal in the different regions.  I see your worksheet does get a solution.  Why do you think it is unrealistic?

Well, cos(m)+I*sin(m)=exp(I*Pi*n) = (-1)^n and cosh(w)+sinh(w)=exp(w) = n^(1/n).
Yes, it's divergent, but the even-numbered partial sums converge to

sum(-(2*n-1)^(1/(2*n-1)) + (2*n)^(1/(2*n)), n=1..infinity)

and the odd-numbered partial sums to that - 1.

I think it has nothing to do with assume.  The difference seems to be the "And".  That is, if any of the conditions in the piecewise is of the form n::something, diff returns unevaluated,
while if the conditions are of the form And(something, somethingelse), it will work.

What do you mean by the derivative with respect to an integer variable?  The derivative as usually defined involves a real or complex variable, although you can also take the derivative with respect to a symbolic variable in a polynomial ring.  But none of these fits what you are doing.  Maple's answer for diff(h1,n) is not correct either.  It should be
piecewise(n::odd, undefined, 3*n^2).

rootlocus(f(s), s,0..30) plots the zeros (in the complex s plane) of 1 + k*f(s) for k from 0 to 30. 

In your case f(s) = 6/(s*(.2*s+1)*(.5*s+1)), and 1+k*f(s) = (s^3+7*s^2+10*s+60*k)/(s*(s+5)*(s+2)).
So there should be three zeros.  Now at k = -14/405+19/810*19^(1/2) ~= .06767787645 the discriminant is 0 and there is a bifurcation, with two real zeros for 0 < k < .06767787645 and two complex conjugate zeros for  k > .06767787645, all coming together at s = -7/3+1/3*19^(1/2) ~= -.880367018. 
I think that, because the plotting is being done using numerical methods, the actual location of the bifurcation is not plotted, and instead you have four points in a diamond-shaped configuration, representing locations of roots just before and just after the bifurcation (in fact, it seems, at k = 15/256 and 5/64).  The other point plotted is actually two points, representing the other root at those two values of k.

This is somewhat tricky, I think, because the left side is 0 at the origin. The curves

y^2 - 2 cos(x) = c

are closed near (0,0).  If I parametrize one of them by theta = arctan(y,x), I find that if u(x,y) <> 0, diff(u,theta) < 0.  Since it is impossible to have this all the way around the closed curve, there are no real solutions in the neighbourhood of the origin other than 0.

Converting the entries to rationals (at Digits = 100, so there's no roundoff), I get the exact value of the determinant:

> with(LinearAlgebra):
   M:= (your matrix):
   Digits:= 100:
   MR:= convert(M,rational):
   Determinant(MR);

-23375023927715521128227661652018007841517045481118503264447374643838762\
207238783235469314216134225716788025159702692516279679196960348785581581\
817776964727838197171425715569931834471184978387901/50000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000

> evalf(%);

-.46750047855431042256455323304036015683034090962237006528894749287677524\
414477566470938628432268451433576050319405385032559358393920697571163163\
635553929455676394342851431139863668942369956775802000000e-5

What you want to modify (in the Classic interface) is the Maple Input style.

Is this what you're looking for?

> U:= u(x,t); 
pde:= diff(U,t$2) = diff(U,x$2);
ic:= u(x,0) = cos(Pi/2*x),D[2](u)(x,0)=0;
bc:= D[1](u)(0,t)=0,u(1,t)=0;
S:= pdsolve(pde,{ic,bc},numeric);
S:-animate(t=0..10);

What are x and y?  What is your code?

If by I you mean the identity matrix rather than Maple's I = sqrt(-1), you might try

> with(LinearAlgebra):
   allvalues(solve(convert(A^2 - 3*A - Delta*IdentityMatrix(2), set),{a,b}));

Note that solve is for finding exact closed-form solutions.  Numerical accuracy should not be a consideration, unless the equations themselves contain floats.  However, if your equations are "quite complex" and non-polynomial in nature, Maple is unlikely to be able to find closed-form solutions. You could try fsolve, which is for finding numerical solutions.

The problem with RealDomain looks like a bug.  Can you give more details, in particular the actual expression final[5]?

On second thought, it is possible to make the number of 1's and -1's equal to 4 with integer linear programming.  For each i, you have two new binary variables a[i] and b[i] with the constraints w[i] = a[i] - b[i], a[i] +b[i] <= 1 (so the possibilities are w[i] = -1 with a[i]=0, b[i]=1, w[i]=0 with a[i]=b[i]=0, w[i]=1 with a[i]=1, b[i]=0), and then
add(a[i] + b[i],i=1 .. NC) = 4.

Your conditions say that each element of A is present once in each row, and once in each column.  This implies that you must have m=n.  Moreover, for any element a of A, the matrix with entries x[i,j] = 1 if a is in m[i,j], 0 if not, is a permutation matrix.  And there is no constraint linking the permutations for different a's.  You could use randperm in the combinat package to construct a permutation for each member of A.  Then the (i,j) element of your Matrix could be the set of all a in A such that the permutation p for a has the property p[i] = j.   You could use select to construct this.

with(StringTools): with(PatternDictionary):
bid:= Create(`builtin`):
words:= [seq](Get(bid,i),i=1..Size(bid)-1):
Anagrams("atc",words);