@Zeineb My suggestion was that you set the values in the pdsolve to match those used in your code. If you've changed the ones in your code, you may have introduced an unrelated problem  

@Zeineb For Crank-Nicolson, it is well known that the oscillations that you describe can occur when a*timestep/spacestep^2 > 1/2, where a is the coefficient of the 2nd derivative in the PDE (I think a=1 in your case; can't check at the moment). So, you need to lower your timestep. See the Wikipedia article on Crank-Nicolson.

@acer Is it safe to put interface(rtablesize= ...) or other interface settings in one's initialization file? The reason that I ask is because I wonder whether the invocation of the initialization file caused by a restart is considered to be in a separate execution group from the restart.

@mmcdara Yes, it should be. The search features of this site are bad, and furthermore, I don't think that Replies are searchable at all.

If you appreciate this Answer, you should give it a vote up.

@Zeineb You asked:

• How can  I use series in Unapply, so that I can display all peicewise functions. If I have two or three for four I can modify by hand but If I have more how can I use many peices in piecewise with unapply

Using the same names that you used in your Question above, that pattern can be continued for an arbitrary number of terms with seq, like this:

unapply(
    piecewise(
        seq(
            eval([x < Xlist[k+1], S[k]], NaturalCoeffs)[], 
            k= 1..numelems(S)-1
        ),
        S[-1]
    ),
    x
);

@Johan159 

Giving a thorough answer requires me to first define for you two of the essential components of function definitions and their invocations. Much of this Answer applies to almost any computer language, although the syntactic minutiae may vary. As always, I use upright boldface (either black or in color) for Maple syntax.

Definition by example of parameter:

• After defining f:= (x,y)-> 2*x^2 + 5*y+3 + 5, the x and y are called the (formal) parameters of f.

Definition by example of argument:

• In these usages: f(3,5), f(x,y), f(a,b)---the a, b, 3, 5, x, y, are called arguments of f (but not necessarily the arguments of f). This is true regardless of whether f has been defined; if it has been defined, regardless of the names of its parameters; and regardless of whether the expression f(...) has been evaluated.

Contrasting parameter and argument:

1. The words parameter and argument are very often incorrectly used interchangeably, and usually this is not a problem because the meaning can be determined from context.
2. Parameter has some other mathematical usages, as I'm sure you're aware, and the two usages do often appear in the same discussion. If it's necessary to be specific, the parameters being discussed here can be called formal parameters. 
3. Argument also has some mathematical usages, but these are so far removed from what were discussing here that it's extremely unlikely to cause confusion. However, if for whatever reason there's a need to be specific, the arguments being discussed here can be called function arguments.
4. It is obvious that a parameter must be a variable and an argument need not be. Maple further restricts parameters to be symbols (which are names without indices) (see help page ?name).
5. This point is highly relevant to your most-recent Reply: Once a variable has been made a parameter of a function^[*1] (or procedure) any usage of that variable's name outside of that function has no effect on its corresponding parameter inside the function. (This is true in almost all computer languages, the exceptions being some old and cruddy ones such as COBOL.) They're completely different variables, stored at different addresses, that only happen to be spelled the same way.

Now let's consider some of your examples:

Example 1, simple evaluation:

You tried:

x:= 0; y:= 5; eval(f(x,y));

You could (and should) simply use 

f(0,5)

In either case, the eval^[*2] does nothing; the evaluation is done for you "for free" in this example. 

Example 2, working with a coordinate slice of a function:

You tried:

y:= 5; plot(f, 0..10);

As stated in point 5 above, using y outside the function has no effect on the parameter y. Since the parameter y has not been given a numeric value, you get the warning message "unable to evaluate the function to numeric values". This message comes from the plot command.

Then you tried:

y:= 5; plot(f(x), x= 0..10);

Although it's not always necessary that every parameter be matched with an argument, if the execution sequence encounters a parameter without a matching argument, you'll get an error message such as "invalid input: f uses a 2nd argument, y, which is missing" (where "f", "2nd", and "y" are replaced with whatever's appropriate for the situation). This message comes not from plot, but from the kernel.

Then you tried:

plot(f, 0..10, 5);

This makes 5 an argument to plot, not an argument to f. The plot command knows that it has no use for a plain number anywhere other than possibly its first argument, so you get the error messge that you did. This error message is explicitly stated as coming from plot.

Then you tried:

plot(f, 0..10, 5..5);

This one fooled plot, because it is possible that a range of two numbers passed as its 3rd argument could be used. Since parameter y still has no value, you get the same message as the first case.

The following ways will work to plot a coordinate slice of f:

1. plot(f(x,5), x= 0..10)
2. plot(f(a,5), a= 0..10)
3. plot(x-> f(x,5), 0..10)

Assuming that a and x do not have assigned values, 1 and 2 do the same thing. And for these to work, indeed the variable used must not have an assigned value (unless, perhaps, that assigned value is another name).

The third way creates a new function (or procedure). Since this procedure has no name, it's called an anonymous procedure. In this case, it makes no difference whether x has an assigned value. When a procedure is created with the arrow, and it has only one parameter and no type declaration, then the parentheses surrounding the parameter are optional.

[*1] The word function also has a meaning specific to Maple. At the moment, it's not necessary that you understand this special meaning. What we're calling "functions" are formally called procedures by Maple, but even Maple itself will sometimes call them "functions" without implying the special meaning of that word that I mentioned in the first sentence of this paragraph.

[*2] There are two essentially different commands named eval. You used it with one argument above. This usage is needed occasionally to force expressions to become "fully evaluated". (If foo differs from eval(foo), then we'd say that foo wasn't fully evaluated.) An example where a beginning user may need this form of eval is 

eval(f);

to show f's definition as a procedure. Only procedures, tables, and modules require this. 

A closely related usage is eval(..., n) where n is a positive integer. This is a rarely used advanced usage, which need not be discussed here.

You may have seen examples of the unrelated other eval command, which is very commonly used, even by beginners, and I think that you may have been trying to emulate this usage. This usage has 2 arguments, the second argument being an equation or a set or list of equations. It's often referred to--incorrectly--as "two-argument eval". Its purpose is to change the values of specific parts of an expression. Some examples are

eval(2*x^3 + 5*y^2, x= 5);
eval(2*x^2 + 5*y^3, [x= 5, y= 3]); #list
eval(2*x^2 + 5*y^3, {x= 5, y= 3}); #set, but does the same thing

Very often a set that can be (and should be) used as the second argument to eval is the return value of another command, such as solve. The right sides of the equations do not need to be numbers; they can be anything. This form of eval is similar to subs. The difference is that eval understands some mathematical subtleties that subs doesn't. The left sides of the equations do not need to be variables (although they are in the vast majority of practical cases).

@fatemeh1090 It can be done, but note that in your plot, x= z, not the Cartesian coordinate x. Anyway,

w:= subs(
    [x= z(r)/sqrt(60), theta= phi],
    0.01503546462*(sin-.0.1328620030*sinh)(-2.365+9.46*x)*cos(6*theta)
        - 0.1
):
(a, R, z):= (2, -8, r-> sqrt(R^2 - (r - a + R)^2)):
domain:= phi= -Pi..Pi, r= 2..3, coords= cylindrical, grid= [50$2]:
plot3d([[w, phi, z(r)], [subs(z= -z, w), phi, -z(r)]], domain);