@acer Expanding on your method, any Unicode "point" in "U+x" form (where x is a string of hex digits) can be converted to a symbol via

`convert/Unicode`:= 
    (U::And(
        string, 
        Not({0,1,2}) &under length,
        identical("U+") &under (index, ()..2),
        satisfies(
            U-> StringTools:-AndMap(StringTools:-IsHexDigit, U[3..])
        )
    ))->
    nprintf(`&#%ld;`, convert(U[3..], ':-decimal', 16))
:
convert("U+03A9", Unicode);
 

@acer Thanks for spotting that, and I have corrected the code above.

@dnaviaux By "style", I only meant with respect to word-processung features such as I mentioned. As shown in my other Answer, adding nonASCII characters to printf output is easy if you're content with fixed character width.

@dnaviaux If that style of report is "fine" with you, then it's trivial to do with printf, and I withdraw my statement about it being "nearly impossible". I was thinking of standard word-processing functionality: word wrapping, page breaks, bold and italics, sub- and superscripts, smoothing the right edges of paragraphs, and supporting non-fixed-width fonts.

I agree with acer: I've changed many of your Replies to Answers and usually given them a Vote Up in the process, all the while thinking "What was he thinking? This is obviously an Answer." But only the one who wrote the Question can select an Answer as Best.

So many of my good Answers don't get Voted Up. It does sadden me. Many users are completely unaware of those buttons. 

In dismantle(a*x + b*y) (as shown by vv above), why are there 5 exponents per term, not 4, and why is the first one 2?

@isabelmacpherson The procedure Kepler below shows how to formulate the system as a vector of four 1st-order ODEs. The rest of the code shows how to perform iterative IVP methods (such as 4th-order Runge-Kutta) using that procedure.

Download RK4.mw

The output is in the downloadable worksheet above. Here's the plaintext form with most output removed:

restart:
Digits:= 15:

#Input ODE IVP system the "normal" way:
xy:= [x,y]: xyt:= xy(t):
Kepler_sys:= {
    (diff(xyt, t$2) =~ -xyt/~add(xyt^~2)^(3/2))[],
    (xy(0) =~ [1,0])[], (D(xy)(0) =~ [0,1])[]
};
                  /  2                                 
                  | d                    x(t)          
   Kepler_sys :=  |---- x(t) = - --------------------, 
                 <    2                         (3/2)  
                  | dt           /    2       2\       
                  |              \x(t)  + y(t) /       
                  \                                    

       2                                                     
      d                    y(t)                              
     ---- y(t) = - --------------------, x(0) = 1, y(0) = 0, 
        2                         (3/2)                      
      dt           /    2       2\                           
                   \x(t)  + y(t) /                           

                             \ 
                             | 
     D(x)(0) = 0, D(y)(0) = 1| 
                              >
                             | 
                             | 
                             / 


#For the sake of comparison, obtain Maple's rk4 solution:
n:= 12:
dsol:= dsolve(
    Kepler_sys, numeric, method= classical[rk4], stepsize= 2*Pi/n
):
#After we get our own solution, we'll compare with Maple's.

#Procedure that evaluates the 1st derivatives. Parameter t is only 
#included for generality; it's not actually used in this ODE system.
Kepler:= proc(t, XY::Vector(4))
local x, y, `x'`, `y'`, d;
    (x, y, `x'`, `y'`):= seq(XY); #Extract input vector.
    d:= -(x^2+y^2)^(3/2);
    #Return the four 1st derivatives (as a vector):
    < `x'`, `y'`, x/d, y/d >
end proc
:
#Procedure that does one step of 4th-order Runge-Kutta:
RK4:= proc(t::complexcons, h::complexcons, Y::Vector, `Y'`::procedure)
local 
    h2:= h/2, t2:= t+h2,
    #All remaining computations use vector arithmetic:
    k1:= `Y'`(t, Y),
    k2:= `Y'`(t2, Y+h2*k1),
    k3:= `Y'`(t2, Y+h2*k2),
    k4:= `Y'`(t+h, Y+h*k3)
;
    Y + h/6*(k1+2*k2+2*k3+k4)
end proc
:   
#This procedure can be used for any IVP method that extrapolates from
#the initial values (which is almost all of them). It defaults to
#using procedure RK4 from above.
IVP:= proc(
    `Y'`::procedure, Y0::Vector[column], 
    T::range(complexcons), n::posint,
    {method::procedure:= RK4}
)
local 
    i, #loop index
    t:= lhs(T), #initial independent variable value
    h:= (rhs(T)-t)/n, 
    Y:= Vector(Y0, datatype= float),
    #matrix for output 
    #(1st row for t-values, 1st column for initial values):
    R:= Matrix((upperbound(Y0)+1, n+1), <t, Y>,  datatype= float)
;
    for i from 2 to n+1 do
        Y:= evalf(method(t, h, Y, `Y'`));
        t:= t+h; 
        R[.., i]:= <t, Y> #Store the ith column of output.
    od;
    R
end proc
: 
Sol:= IVP(Kepler, <1, 0, 0, 1>, 0..2*Pi, n);

#Index [.., -1] is the last column of a matrix.
Sol[.., -1];
                 [     6.28318530717959      ]
                 [     0.991508429966570     ]
                 [    0.0430479365430936     ]
                 [    -0.0436637679707416    ]
                 [     1.00287851473040      ]

#Compare with dsolve's result:
eval(<t, x(t), y(t), diff(x(t),t), diff(y(t),t)>, dsol(2*Pi));
                 [     6.28318530717958      ]
                 [     0.991508429966570     ]
                 [    0.0430479365430934     ]
                 [    -0.0436637679707406    ]
                 [     1.00287851473040      ]

#You should redo with much larger n.

#The following will do a formatted print of every mth t, x, and y value,
#just like you were doing:
m:= 2:
printf("%10.4f\n", Sol[..3, [seq(m..upperbound(Sol)[2], m), -1]]^+);
    0.5236     0.8660     0.4996
    1.5708    -0.0003     0.9981
    2.6180    -0.8646     0.4935
    3.6652    -0.8531    -0.5088
    4.7124     0.0258    -0.9922
    5.7596     0.8777    -0.4635
    6.2832     0.9915     0.0430

@AHSAN I agree with acer 100%. No useful purpose would be served by explicit symbolic representation of the parameterized roots.

Yes, I do see the issue of quotes as a major drawback to this method that I didn't originally anticipate. Almost any Maple code of significant length will contain pairs of both single and double quotes, although quotes inside quotes are a rarity.

@mmcdara What you showed regarding Complex(a, b) is interesting in its own right, but it sn't the point that I was trying to make in my last paragraph. I was referring to Complex(a, b) where a and b are manifest real constants. For example:

ToInert(Complex(3,7));
             _Inert_COMPLEX(_Inert_INTPOS(3), _Inert_INTPOS(7))

So, we see that the result has nothing to do with either functions or names.

You said "running the code"; however, it's not clear to me that that's what you really mean. Are you instead referring to parsing errors such as missing semicolons and parentheses?

The techniques that you mentioned don't apply to basic syntax errors---such as missing semicolons or parentheses---that prevent the code from even being "compiled" or parsed or evaluated. I think that those are the type of errors that the OP is referring to, although I'm not sure, and I've asked for clarification.

Nonetheless, all that you wrote is valuble information to know once one has actual runtime errors.

The (*  *) style comments work fine in Maple 2020. However, they can't span multiple execution groups in a worksheet, which is what I suspect you're trying to do.

@emendes Interesting.... I suspect a bug in Threads. Could you send me (either email or post here) the file that contains the lists of sets A and B for this? I don't think that there's any hope of finishing the computation on my small computer; I just want to see if I can get a strange error. I hope that you can send me a file that contains only two lists of sets that I can read in with the read command.

@aoakindele Simply change your dsolve command to

sol1[k]:= dsolve(eval(odeSys union bcs1, pList), numeric, method= bvp[midrich]);