@nm
solve evaluates its arguments before proceeding as is the general rule in Maple.
So try these:

restart;
x>a and x<b and x<>c;
x>a and x<>c;

restart;
x>-infinity and x<infinity and x<>1/2;

The output discards the '<>' part in both cases.  Thus solve never sees x<>1/2.

In the help page for solve/details we find under Parameters that the first argument can be an
"expression, equation, inequality, or procedure; or set or list of expressions, equations, or inequalities".
Admittedly the term 'expression' suggests something vague.

## Whether the 'generic' handling of the case 'x>a and x<>c' is a bug or not I shall find out when I report it as such by submitting anSCR.

@samen Yes, indeed.
With a few changes here is my version:
 

restart;  #Using the default setting of Digits (10, why lower?)
# Parameter values:
a := 0; b := 1; N := 9; 
h := (b-a)/(N+1); phi := .5; 
K := 10^(-6); mu := 1.67; 
alpha := K/(phi*mu); lambda := alpha*k/h^2;
# Initial conditions
for i from 0 to N do 
  u[i, 0] := h*i+1 
end do:
# Boundary conditions
for j from 0 to N+1 do 
  u[0, j] := .1; 
  u[N+1, j] := .5 
end do:
#Discretization Scheme
for i to N do 
  for j from 0 to N do 
    eq[i, j] := lambda*u[i-1, j]+(2-2*lambda)*u[i, j]+lambda*u[i+1, j] = -lambda*u[i-1, j+1]+(2+2*lambda)*u[i, j+1]-lambda*u[i+1, j+1] 
  end do 
end do:
sys := ([seq])(seq(eq[i, j], j = 0 .. N), i = 1 .. N):
nops(sys); # 90
vars:=indets(sys) minus {k}:
nops(vars); # 90
nn := Matrix(N+1, N+1,(i, j)-> u[i-1, j-1]):
##
p:=proc(kk) local u_res,A; 
  u_res:=solve(eval(sys,k=kk),vars);
  A:=eval(nn,u_res);
  plots:-matrixplot(A)
end proc;
##
p(10000); # Test
plots:-animate(p,[k],k=0..10000,orientation=[145,70,15]);

@samen You could use matrixplot from the plots package:
 

u_res := solve(sys);
nn := Matrix(N+1, N+1,(i, j)-> u[i-1, j-1]);
A:=eval(nn,u_res);
plots:-matrixplot(A);

### I understand that you also want to see the dependence on k.
## Note: lambda is the only actual parameter in the equations. Since h is the space step I suppose that k may be the time step. If so it couldn't be zero. But you also have a K, phi, and mu in the definition of lambda. So varying any one of these in an appropriate interval will produce the same animation as we get by varying k.
##
Here is a way using animation in k.
Delete the assignment to k. Begin with a restart at the top.
Do the boundary and initial value assigments as before.
Define the equations eq[i,j] as before. They will now contain the parameter k (unassigned).
Then do:

sys := ([seq])(seq(eq[i, j], j = 0 .. N), i = 1 .. N):
nops(sys);
vars:=indets(sys) minus {k}:
nn := Matrix(N+1, N+1,(i, j)-> u[i-1, j-1]):
##
p:=proc(kk) local u_res,A;
  u_res:=solve(eval(sys,k=kk),vars);
  A:=eval(nn,u_res);
  plots:-matrixplot(A)
end proc;
## Testing p for k=10000:
p(10000);
## Animating the plot for k=0..10000:
plots:-animate(p,[k],k=0..10000);

@basha 666 I noticed two more problems. The code
 

Rf:=diff(f[m-1](x),x,x,x)+2*alpha*r*sum*(f[m-1-n](x)*diff(f[n](x),x),n=0..m-1)
+(4-Ha)*(alpha)^2*diff(f[m-1](x),x);

has a `*` after 'sum'. That certainly shouldn't be there.
Secondly, in the code:
 

h*H(diff(f[m-1](x),x,x,x))

you have H being treated like a function (i.e. procedure). Since H:=1 is set at start this means that
h*H(diff(f[m-1](x),x,x,x)) evaluates to h.
If you meant multiplication i.e. wanted h*H*(diff(f[m-1](x),x,x,x)) then the expression would evaluate to h*diff(f[m-1](x),x,x,x).
So there is a difference.
These problems may be due to your use of 2D input (do you use that?). But as it appears in your code it makes quite a difference in the more reliable 1D input (a.k.a. as Maple input).

If interface(ansi=true) is set then showstat prints real bad in Windows 10 in Standard Worksheet Maple. In cmaple it looks fine:

restart;
interface(ansi=true);
showstat(apply);
interface(ansi=false);
showstat(apply);

The first showstat shows this:

apply := proc(p)
local P;
   1   eval(subs(P = p,P(args[2 .. -1])))
end proc

I chose the really short procedure apply for convenience.

So in Windows 10 the maple.ini file will have to distinguish between handling Command-line and Standard Worksheet Maple.  Here sscanf can work on interface('version') where the interface used can be found.

@nm The help page referred to by Tom Leslie is ?copy, not ?Record.

@Joe Riel Thank you for sharing this code. I'll be using sscanf from now on.

@tomleslie It is wrong in general. Consider e.g. this:
 

restart;
expr1:=sqrt(x^2+y^2)/x;
expr2:=sqrt(collect(expr1^2,x));
eval(expr1,{y=1,x=-2});
eval(expr2,{y=1,x=-2});

@acer Paulina Chin wrote on December 3 2010:

" The problem with ranges having more than 3 dots that acer and Joe had commented on is indeed a bug and it's already in our database. "

I suppose that that bug is enjoying a long and cosy life in that database.

About putting interface(ansi=true) in my maple.ini file:
It works fine in Maple 2018.1, but in earlier versions it has the effect of making output look real weird.
Thus I now use a construction I already had in place anyway to use for setting ansi=true only in versions 2018.1 and up:
 

proc() local VERSION,ns,N;
VERSION:=convert(kernelopts(version),string);
ns:=StringTools:-Search(",",VERSION); 
N:=parse(VERSION[7..ns-1]);
if N>=2018.1 then interface(ansi=true) end if;
NULL
end proc();

It seems to me to be a bug in the 2D parser. Spaces are not necessary in 1D.
If you want seq(-3...3) to mean seq(-3..0.3) and don't want to write the zero you can insert a space: seq(-3.. .3).
Similarly this returns floats: seq(-3. .. .3);

@stefanv Thank you. I had success with "Standardindstillinger" rather than "Egenskaber" (in English presumably "Preferences" and "Properties", respectively). With "Properties" the font size changed right away, but an error message showed that a shortcut couldn't be established. "Preferences" worked as you described.

I haven't been using the command line interface except for trying it as I just did again.
The first that strikes me in Windows 10 is the ridiculously small text size used. That is enough for me to quit cmaple.
I suspect that it can be made bigger, but I don't know how.

You say that highlighting is the default in Windows. Not so on my version. The command
interface(ansi);
returns false.

@nm You can use combine followed by collect:
 

restart;
f := x->3*cos(42*x*Pi);
g :=(x,t)->exp(3*t)*cos(17*x*Pi);
pde := diff(u(x, t), t)= (diff(u(x, t), x, x)) + g(x, t);
bc := (D[1](u))(0, t) = 0, (D[1](u))(1, t) = 0;
ic:= u(x, 0) = f(x);
sol:=pdsolve([pde, ic, bc]);
##
combine(sol);
collect(%,exp,factor);

@panke The code works just the same in Maple 12 (I don't have Maple 13 available).
The line

SOL(eventfired=[1]);

can be left out if you wish. It just shows the same result basically.