@mmcdara 

I just extensively updated the code and examples in my most-recent Reply, including the section that you just asked about. Please read that before continuing here.

Your question is very good, but it deals with one of the most-arcane aspects of Maple syntax: the overloading of infix operators. So, I'm sure that after reading this, you'll have several further questions. Please feel free to ask those questions.

Like most Maple operators, the elementwise operator ~ is overloadable. That means that its meaning can be changed for certain cases without changing its overall meaning. The filtering of the "certain cases" is determined by passing or failing the type checks of the procedure headers of a list of procedures in an overload command. (Operators acting on objects can also be overloaded, but by a completely different mechanism that isn't used here.)​​​​​​

Akin to most infix operators, the procedure that controls ~ is named `~`. However, it's a bit more arcane than most other operators because the procedure itself is invoked with an index. Specifically, the code f~(a, b, c) invokes `~`[f](a, b, c). That's why the first procedure in my overload uses op(procname), which extracts the index from an indexed procedure invocation. Some further arcanity is needed if f is builtin, in which case I replace f with f@(x-> args) in the overload.

@mmcdara You wrote:

• It could be easy to transform Kitonum's output or mine into a matrix by doing `<|>`(op(%))​​​​​​

That's not automatic: It only works on matrices (as opposed to higher-dimensional arrays) and only when the operation is columnwise.

• Your procedure is indeed an interesting to apply a function f : U-->V where U and V are vectors of the same length.

The generalized procedure below allows for them to be not of the same length; indeed, the Vs can be anything at all, and they need not even be of the same type or size. And the Us need not be vectors; they can be any "slice" (or subarray) of an input Array, of any number of dimensions.

• Could you adjust in such a way it also supports simpler transformations, for instance to mimic abs~(v), or cos~(v) , in order to make it a more versatile tool?

I don't know why you'd want that given that abs~(v) and cos~(v) already work as is; nonetheless, the code below handles that case. Specifically, if the set in the argument dims={...} contains all the dimensions, then the result of DimMap (whether it's invoked directly or via ~) should mimic ordinary elementwise operation.

restart:

#Decide on the rtable subtype of the result: Vector[row], Vector[column],
#Matrix, or Array:
 
Subtype:= proc(A::rtable, dims::set(posint))
    `if`(
        A::Vector,
        rtable_options(A, ':-subtype'),
        `if`(
            A::Matrix,
            `if`(
                nops(dims)=1,
                Vector[`if`(dims={1}, ':-column', ':-row')],
                Matrix
            ),
            Array
        )
    )
end proc
:
#Main procedure:
DimMap:= proc(f, A::rtable, Dims::{identical(dims)=set(posint)})
local
    j, _j, 
    rest:= _rest, 
    i:= 0, inc:= proc() i:= i+1 end proc,
    d:= rhs(Dims), 
    D:= [rtable_dims](A),
    _J:= [seq](`if`(j in  d, _j[inc()], D[j]), j= 1..nops(D)),  
    F:= rtable(
        D[[d[]]][], 
        ()-> f(A[eval(_J, _j= [args])[]], rest),
        ':-subtype'= Subtype(A, d)           
    )               
;
    if F::'rtable'(rtable) and (nops@{entries}@rtable_dims~)(F) = 1 then
        rtable(evalindets(F, rtable, convert, list, ':-nested'))      
    else
        F
    fi
end proc
:

#Overload elementwise operator ~ to work with DimMap:

unprotect(`~`, `~orig`):
`~orig`:= eval(`~`):
`~`:= overload([
    proc(A::rtable, Dims::{identical(dims)=set(posint)}) 
    option overload;
    local f:= op(procname);
        DimMap(`if`(f::builtin, f@(x-> args), f), args)
    end proc,
    `~orig`
]):
protect(`~`, `~orig`):

#Examples:
Sx:=1/sqrt(2)*Matrix([[0,1,0],[1,0,1],[0,1,0]]):
lam,v:=LinearAlgebra:-Eigenvectors(Sx);

                              [-1]  [   1     -1    1   ]
                              [  ]  [                   ]
                              [ 0]  [  (1/2)       (1/2)]
                    lam, v := [  ], [-2       0   2     ]
                              [ 1]  [                   ]
                                    [   1     1     1   ]

#Normalize columnwise in 2-norm. (In the following 4 examples, the final 2 could be
#replaced by 'Euclidean' and the results would be the same.)
LinearAlgebra:-Normalize~(v, dims={2}, 2);

                     [    1         1  (1/2)     1    ]
                     [    -       - - 2          -    ]
                     [    2         2            2    ]
                     [                                ]
                     [  1  (1/2)              1  (1/2)]
                     [- - 2           0       - 2     ]
                     [  2                     2       ]
                     [                                ]
                     [    1        1  (1/2)      1    ]
                     [    -        - 2           -    ]
                     [    2        2             2    ]

#Normalize rowwise in 2-norm:
LinearAlgebra:-Normalize~(v, dims={1}, 2); 

                     [ 1  (1/2)     1  (1/2)  1  (1/2)]
                     [ - 3        - - 3       - 3     ]
                     [ 3            3         3       ]
                     [                                ]
                     [  1  (1/2)              1  (1/2)]
                     [- - 2           0       - 2     ]
                     [  2                     2       ]
                     [                                ]
                     [ 1  (1/2)    1  (1/2)   1  (1/2)]
                     [ - 3         - 3        - 3     ]
                     [ 3           3          3       ]

#Compute 2-norms columnwise:
LinearAlgebra:-Norm~(v, dims={2}, 2);

                               [    (1/2)   ]
                               [2, 2     , 2]

#Compute 2-norms rowwise:
LinearAlgebra:-Norm~(v, dims={1}, 2);

                                  [ (1/2)]
                                  [3     ]
                                  [      ]
                                  [  2   ]
                                  [      ]
                                  [ (1/2)]
                                  [3     ]

#Compute abs elementwise:
abs~(v, dims={1,2}); 

                             [  1     1    1   ]
                             [                 ]
                             [ (1/2)      (1/2)]
                             [2       0  2     ]
                             [                 ]
                             [  1     1    1   ]

#Multi-sub-dimensional examples: A is 2x2x2. Indexing along the
#2nd dimension, we multiply all entries in the 1st and 3rd dimensions.
#So, the result is a 2-element vector. Then we multiply between those
#same dimensions, producing a 2x2 result.

A:= Array((1..2)$3, rand(1..8)): 
'A'[.., 1, ..] = A[.., 1, ..], 'A'[.., 2, ..] = A[.., 2, ..];

`Product of slices in the 2nd dimension` =
    (`*`@op@op~@[entries])~(A, dims={2})
;

`Product between those same slices` =
    (`*`@op@op~@[entries])~(A, dims={1,3})
;

                               [3  4]                             [4  6]
    A[() .. (), 1, () .. ()] = [    ], A[() .. (), 2, () .. ()] = [    ]
                               [5  3]                             [8  1]
             Product of slices in the 2nd dimension = [180, 192]
                                                    [12  24]
                Product between those same slices = [      ]
                                                    [40   3]

#In recent versions of Maple (that support one-argument mul), the arcane operator
#(`*`@op@op~@[entries]) could be replaced by simply mul:

mul~(A, dims={2}), mul~(A, dims={1,3});
#

Several points:

1. What you call "factors" of a natural number are properly called divisors. However, I think that your intended meaning is clear. If the divisors are prime, then it's okay to call them factors, but it'd be better to be specific and call them prime factors.

2. On page 34, you say "[The] Euler totient function is the number of factors of a natural number." That is completely wrong. Indeed, it's almost the opposite of that: The Euler totient of n is the count of natural numbers a less than n such that gcd(a, n) = 1. So, it's a count of numbers that share no divisors (other than 1) with n.

The count of the positive divisors of n is often called sigma₀(n).

3. Although I read your PDF document quickly, I think that several pages are duplicated. Some other pages seem to be out of logical sequence.

4. Any elementary presentation on prime-generating polynomials should include a proof (since the proof is elementary) that no polynomial can produce prime output for all natural-number input.

@HebaSami 

This is a direct Answer to your titular question "How do I use numerical solutions of [ODEs] in other equations ...?": The process that I showed above will work to produce plots of any expressions created from the dependent variables of an IVP as long as those expressions don't contain integrals and don't contain derivatives of the same or higher order than the highest-order derivatives in the ODEs. With some small modifications, any of the following can be handled:

• BVPs,
• higher-order derivatives,
• numeric integration of expressions of the dependent variables,
• end results other than plots that require numeric evaluation of the expressions.

@tomleslie Thank you, Tom. It seems that I've become habituated to the new seq syntax. By the way, you may realize that the only reason that I put those expressions in vectors is to improve their prettyprinted display. 

@HebaSami The error message suggests to me that you didn't transcribe the angle brackets <  > that I put in the command that defines Extra. If you can't fix it, upload an executed worksheet containing the error message. Use the green up-arrow on the editor toolbar. Ignore any error message that the uploader gives. 

@MuYi I'm curious whether Maple can access your swap. Issue the command

kernelopts(datalimit);

The value returned is the number 2^10 = 1024-byte blocks that Maple can access. You may need to increase this. If so, let us know.

@ The process that you've elaborately described as "goal seeking" is commonly called "root finding". It is one of the fundamental problems of numerical analysis, and one of the oldest problems in mathematics. There are at least tens of algorithms and thousands of computer programs devoted to it. If you're interested in studying numerical analysis, then we can discuss some of those algorithms and write some procedures. But if you're only interested in a practical solution to your problem, then you're wasting your time writing your own algorithms. Although your function G may seem complicated to you, it is in fact a very simple root-finding problem. So, if the existing algorithms can't find the root, it likely can't be found at all.

The main Maple command for root finding is fsolve, but there are also several other commands. The version of G that you've posted immediately above can be easily solved by fsolve:

restart:
Digits:= 15: #optional
G:= eval[recurse](
    R*T*rho*ln((s + 1)*(1 - rho*(1-s))/((1-s)*(1 - rho*(s+1)))) 
        - 4*s*rho*(1 - rho)*((1 - rho)*A + rho*B),
    [R= 8.314, A= 3.511*R, B= 4.389*R, T= 1.0, rho= (T/3.135)^(3/2)]
);
r1:= fsolve(G, s= 0..1);
                            r1 := 0.

#I'll assume that you're not particularly interested in that root,
#so use the 'avoid' option to find another. This is the general way
#to find multiple roots with fsolve.

r2:= fsolve(G, s= 0..1, avoid= {s=r1});
                    r2 := 0.999981403075794

#It can be easily seen from a plot that these are the only two roots.

Note that the initial value of s that you checked in your procedure was 0.995, which is already significantly less than the true root.

If you have some other version of G for which fsolve is not working for you, please post it.

@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);