@ I don't see any attached file. Please try again. 

The following symbols are used in your PDF without the slightest definition, not even an equation; nor are they used in your worksheet, so that can't be used to figure out their meanings either: x, y, a, v__f, w, C, T, N.

"In the ... absence of lambda" makes no sense me, unless perhaps by "absence" you mean lambda=0.

I can see that the command Embed exists, but it is undocumented. What is it supposed to do?

@Carl Love Here's a combined plot of the entries computed by the "novel rule" and by the "repeat rule". It shows an interesting pattern.

A:= proc(N::And(posint, Not({1,2})))
local 
    n, m, p, U:= table([HFloat(1)= 1]), 
    B:= rtable(1..N, [false$2], datatype= truefalse), 
    A:= rtable(1..N, [1,0], datatype= float[8])
;
    for n from 2 to N-1 do
        p:= A[n];
        m:= U[p];
        B[n+1]:= type(m, posint); 
        A[n+1]:= `if`(B[n+1], n-m, min(abs~(p -~ A[..n-1])));
        U[p]:= n       
    od;
    ([seq](trunc(a), a= A), B)
end proc
:
N:= 2^13:
(S,B):= CodeTools:-Usage(A(N)):
memory used=152.32MiB, alloc change=27.49MiB, 
cpu time=469.00ms, real time=491.00ms, gc time=0ns

#Separate the indices into 2 lists: "repeat rule" and "novel rule":
(i1,i2):= selectremove(k-> B[k], [$1..nops(S)]):
evalf(nops(i1)/N);
                                0.7999267578

plot(
    [`[]`~(i1, S[i1]), `[]`~(i2, S[i2])],
    style= point, legend= [Repeat, Novel], color= [red, blue],
    symbolsize= [1,9], axes= boxed, symbol= [point, solidcircle]
);

@David Sycamore The sequence output by the procedure is sufficient for either of those two purposes. There is nothing that needs to be "adapted".

To plot: Using the sequence S already returned, do plots:-listplot(S, ...) where ... are options exactly as shown in the other Answers.

To check that inequality:

remove(n-> S[n]+S[n+1] <= n, [$1..nops(S)-1]);
               [521, 1792, 2377, 3670, 4332, 5767, 6786, 7442]

The returned list is the values of n for which the inequality is false.

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