I'm glad it brought that to your mind, Preben, I'd forgotten about it completely.  For those interested, here it is (from comp.soft-sys.math.maple in November 2001:
<http://groups.google.com/group/comp.soft-sys.math.maple/msg/66f768adf8827e69?hl=en>):

>   applysome:= proc(fn, S::set, expr) 
    local yes, no,i,n, f0; 
    n:= nops(expr); 
    f0:= op(0,expr); 
    yes:= NULL; no:= NULL; 
    for i from 1 to n do 
      if member(i,S) then 
        yes:= yes,op(i,expr) 
      else no:= no, op(i,expr) 
      fi 
     od; 
    f0(no,fn(f0(yes))) 
    end;

Nowadays I wouldn't use a for loop:

> applysome:= proc(fn, S::set, expr)
    local yes, no, f0;
    f0:= op(0, expr);
    yes:= seq(op(i,expr), i=S);
    no:=  seq(op(i,expr), i={$1..nops(expr)} minus S);
    f0(no, fn(f0(yes)))
    end proc;

Thus applysome(fn, S, e1*e2*...*en) will return the product of those ej
for j not in S and fn of the product of those ej for j in S.  Thus

> applysome(expand, {2,3}, (x-1)*(x-2)*(x-3)*(x-4)); 

(x-1)*(x-4)*(x^2-5*x+6)

For another example:

> applysome(`*`@op, {1,3,5}, a1+a2+a3+a4+a5); 
 

                       a2 + a4 + a1 a3 a5

I'm glad it brought that to your mind, Preben, I'd forgotten about it completely.  For those interested, here it is (from comp.soft-sys.math.maple in November 2001:
<http://groups.google.com/group/comp.soft-sys.math.maple/msg/66f768adf8827e69?hl=en>):

>   applysome:= proc(fn, S::set, expr) 
    local yes, no,i,n, f0; 
    n:= nops(expr); 
    f0:= op(0,expr); 
    yes:= NULL; no:= NULL; 
    for i from 1 to n do 
      if member(i,S) then 
        yes:= yes,op(i,expr) 
      else no:= no, op(i,expr) 
      fi 
     od; 
    f0(no,fn(f0(yes))) 
    end;

Nowadays I wouldn't use a for loop:

> applysome:= proc(fn, S::set, expr)
    local yes, no, f0;
    f0:= op(0, expr);
    yes:= seq(op(i,expr), i=S);
    no:=  seq(op(i,expr), i={$1..nops(expr)} minus S);
    f0(no, fn(f0(yes)))
    end proc;

Thus applysome(fn, S, e1*e2*...*en) will return the product of those ej
for j not in S and fn of the product of those ej for j in S.  Thus

> applysome(expand, {2,3}, (x-1)*(x-2)*(x-3)*(x-4)); 

(x-1)*(x-4)*(x^2-5*x+6)

For another example:

> applysome(`*`@op, {1,3,5}, a1+a2+a3+a4+a5); 
 

                       a2 + a4 + a1 a3 a5

@jean-jacques : hirnyk's code works in 1D input.  You appear to have typed sea instead of seq.  Try copying and pasting rather than typing it yourself.

@jean-jacques : hirnyk's code works in 1D input.  You appear to have typed sea instead of seq.  Try copying and pasting rather than typing it yourself.

The constant is still the same constant.  Values of u or u[t] (at the mapped values of X, of course) are the same.  The only things that must be scaled are derivatives with respect to x.

The constant is still the same constant.  Values of u or u[t] (at the mapped values of X, of course) are the same.  The only things that must be scaled are derivatives with respect to x.

I don't see why not.

I don't see why not.

Here is an example.

Xpde.mw

Here is an example.

Xpde.mw

The curvature of the graph of y = f(x) is abs(diff(f(x),x$2))/(1+diff(f(x),x)^2)^(3/2).

>

The curvature of the graph of y = f(x) is abs(diff(f(x),x$2))/(1+diff(f(x),x)^2)^(3/2).

>

What varies in your PDE is the thermal diffusivity.  The thermal conductivity does not appear in the PDE.  If you just solve the PDE with the piecewise thermal diffusivity, you are in effect assuming that the thermal conductivity is constant, so that the change in thermal diffusivity is due to a change in the heat capacity.  That leads to the results you obtained.

For example, suppose for a[i] <= x <= a[i+1] you have thermal diffusivity r[i] and thermal conductivity k[i]  We can deal with this by a piecewise linear change of the spatial variable.  If X is the new spatial variable, you want Delta X = Delta x/k[i] in region number i, and
the PDE there is u[t] = r[i]/k[i]^2 * u[XX]. 

What varies in your PDE is the thermal diffusivity.  The thermal conductivity does not appear in the PDE.  If you just solve the PDE with the piecewise thermal diffusivity, you are in effect assuming that the thermal conductivity is constant, so that the change in thermal diffusivity is due to a change in the heat capacity.  That leads to the results you obtained.

For example, suppose for a[i] <= x <= a[i+1] you have thermal diffusivity r[i] and thermal conductivity k[i]  We can deal with this by a piecewise linear change of the spatial variable.  If X is the new spatial variable, you want Delta X = Delta x/k[i] in region number i, and
the PDE there is u[t] = r[i]/k[i]^2 * u[XX]. 

The inverse of the "exact" version of your Matrix has some huge entries, on the order of 10^35, so it's not surprising that you get numerical difficulties.   This sort of thing can happen with ill-conditioned matrices. You didn't tell us what problem this came from (in particular, you only gave us one 6 x 6 matrix, so we have no idea what is happening with N>=7.  I wonder if what you have is an approximation of a "real" matrix that is singular.