Dear all, as a former aerospace engineer working in the field of astrodynamics and orbital mechanics I had a great need to work with Julian Dates and converting them back to a Gregorian Calendar dates.  The source of this algorithm I use is from the book by David A. Vallado "Fundamentals of Astrodynamics and Applications" his 3rd edition however he is up to his 5th edition and the algorithm is still there.  I had to make a change to his algorithm to deal with a leap years and the months of the year with only 30 days, eg., April, June, and November.  I know this is not as elegant as acer's and others but I have used it for many years and I am satisfied it works!  Note that in Maple, if you use the JulianDayNumber from the Calendar package YOU NEED TO SUBTRACT 0.5 from the number it produces!!! This is because astronomers don't count days from high noon but from midnight.  They are lazy and that is their convention for the start of a new day! If you want to check this algorithm, look at any edition of the "Astronomical Almanac" published by the U.S. Naval Observatory and they list the JD for each day of the year to which that almanac is published for.  It is published yearly.   

JDToGregorianDate(0.24498773458762e7, SecDecimals = 6, TABOUT = true);
JD (UTC)
Year
UTC
Month
UTC
Day
UTC
Hours
UTC
Minutes
UTC
Seconds
UTC
Day Of Week
UTC
2449877.3458762
                             "1995"

                              "6"

                              "8"

                              "20"

                              "18"

                            "14.400"

                           "Thursday"

Download JDToGregorianDateConverter.mw 

Hello! I once use to controls system analysis and when we needed to represent a delay, we used Pade approximations for a delay.  If you want to see this, I will search my Maple files for such an example.

I have done some more work on my Mars flyby trajectory Maple worksheet.  

Thank you again for your help on the piecewise function problem I was having.

@SUN ZHUOYU You can adjust the fonts on the axes using, for example, axesfont=[Helvetica,bold,18].  Also, again for example,the label's font using, for example, labelfont=[Helvetica,bold,24] inside your display command.

q1_(1).mw

Take a look at this and see if it is what you expect.  Perhaps you can "eye-ball" the maximum value for a particular q and p.

plot(rhs(solution), t = 0 .. 600, view = [default, -600000 .. 0], size = [1200, 1200], gridlines = true, color = red, thickness = 5, background = "Azure", axesfont = [Helvetica, bold, 24])