@schapplm There is a procedure sqrt in Maple. You can see the code by doing
 

showstat(sqrt);

In fact it appears that at least since Maple 15 sqrt is actually a module that has ModuleApply as a local procedure. This means that the module sqrt can be used as a procedure. That procedure is the procedure ModuleApply.
Details:
 

restart;
showstat(sqrt);
kernelopts(opaquemodules=false);
showstat(sqrt:-ModuleApply); # same procedure as the one above
op(3,eval(sqrt:-ModuleApply)); # option builtin, but can be seen
op(2,eval(sqrt)); # locals include 3 procedures and one large integer.
pr:=sqrt:-Primes; # The product of the first 35 primes:
ifactor(pr);
nops(%);
ithprime(35);
showstat(sqrt:-Product); 
showstat(sqrt:-Nested); 

@tomleslie Yes, just by replacing fsolve by solve (and using sigma1=1/2) makes sol = NULL.

About the boundary conditions: You have
D(f)(0)=1, D[1,1](f)(0), D[1,1,1](f)(6)=2
in your call to dsolve, but in the discrete version we find
 

bcEqs:=[subs(k=0,dfdef(h))=rhs(bc1),subs(k=0,d2fde2f(h))=rhs(bc2),
subs(k=n-1,d2fde2b(h))=rhs(bc3)];

which are discretized versions of D(f)(0) = rhs(bc1), (D@@2)(f)(0)=rhs(bc2), (D@@2)(f)(6) = rhs(bc3).
The notations d2fde2f and d2fde2b refer to forward and backward discretized second derivatives, respectively.

On my Windows 10 computer I get that sound when using the arrow keys to move from one input line to the next but only when there are no more input lines in the direction of the arrow.
That doesn't annoy me and is nothing new. I have the sound as early as Maple 12 (Standard Worksheet Interface), but not in Maple 8 (Classical Worksheet).
Could you give an example?

## I can produce the same sound by using the Delete key or the back delete key (or whatever it's called) when there is nothing to delete.
 

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