@cky1946 You can find the information on the help page ?plot,options. (There's a link to this from ?plots,matrixplot.) This is one of my most frequently reread help pages.

@tomleslie Thank you, Tom. I think that most of the more-abstract plotting commands should wrap their output with a display command such as plots:-display(..., overrideoption, _rest) so that any arguments not matching declared keyword parameters are just passed through to display.

@tomleslie Tom: Your confusion arises because although the kernel only "sees" 1-byte characters, the Maple user interfaces (including the plaintext command-line interface) do fully support multi-byte-encoded characters.

Vote up for an impressive solution. 

If you wish, your cat command can be replaced by sprintf("%a[%d:%d]", s, m-1, n). I find that easier to read, although I do realize that what's easier for me to read often doesn't correspond to what others find easier to read.

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