@tomleslie I'm using 64-bit Windows 10.

@farah adanan In the paper you link to the odes (13), (14), and (15) together with the boundary conditions (16) really will be a boundary value problem, thus you will not be using Runge-Kutta-Fehlberg.
Infinity in (16) must be replaced by a finite number of an appropriate size (not huge).
Mapleprimes has had quite a lot of problems in that ball park.
You should try to write the equations using Maple syntax. Then we could take it from there.
Boundary value problems are often difficult to solve.

@tomleslie For the same commands I get:

CUDA:-ComputeLevel(id=0);
                               5.
CUDA:-HasDoubleSupport(id=0);
                              true
CUDA:-Properties();
                              [table(["Total Global Memory" = 2147483648, "Minor" = 0, "Name" = "GeForce GTX 860M", "Texture Alignment" = 512, "MultiProcessor Count" = 5, "Device Overlap" = 1, "Clock Rate" = 1019500, "Max Threads Per Block" = 1024, "Max Grid Size" = [2147483647, 65535, 65535], "Max Threads Dimensions" = [1024, 1024, 64], "Major" = 5, "Memory Pitch" = 2147483647, "Resisters Per Block" = 65536, "Total Constant Memory" = 65536, "Kernel Exec Timeout Enabled" = true, "Shared Memory Per Block" = 49152, "Warp Size" = 32, "ID" = 0])]

@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.