What do the regions labelled  TM1, TM2, TM3 are aimed to represent?

Did you mean this: Regions_where_each_TMi_is_the_largest.mw

Download One_way.mw

Something_like_that_maybe.mw

(Shame there is no example in the corresponding help page about this way to use dualaxisplot)

>  >  >  >  >  (1) >  (2) >  (3) >  >  (4) >  (5) >  (6) >  >  dualaxisplot(   display(S1, S2, S3)   ,   display(S12, S22, S32) )

Download All_plots_Combined_sand15.mw

One way is to give rho a numeric value BEFORE solving the ode system, the other one (which I use in the attached file) is to use the dsolve/numeric option 'parameters'.
 

Before doing this I advice you to rewrite the system you want to solve, which is actually a DAE (Differential Algebraic System), as an ODE system (Maple can solve DAE systems but I can't see any reasons for doing so here unless adding unnecessary complexity).
 

Once all this done dsolve is then capable to compute a solution. 
Nevertheless you will see in the attached file that an arbitrary choice of the rho value (I took rho = 1), leads to a singularity at time 0.57: because I have no idea of the values rho might have

(rho is ought to represent the density of Mars atmosphere, something about 20g/m³, but I don't know what are the units you use) 

it's hard to say more about this dicontinuity (real or simple artefact due to unrealistic rho value?).
Nevertheless a quick analysis is provided at the end of this answer.

SimpleMarsEntryAndAeroBrakingModel_sand15.mw



Note that other choices enable computing the solution for arbitrary much larger times, for instance

ans(parameters=[0.1])
                          [rho = 0.1]

display(
  Matrix(2, 3, [seq(odeplot(ans, [t, u], t=0..100, title=typeset(u)), u in unknowns)])
)

Nevertheless I suppose you want to end the simulation when the probe hits the ground, don't you?  
If it is so, let landing_t the time when the probe reaches Mars surface. Given its average radius is 3389.5 Km, landing_t verifies  r(landing_t) = 3389.5. The easiest way to stop the simulation when this latter condition is fullfiled uses the dsolve/numeric 'events' option:

landing_r := 3389.5; # kilometers

ans := dsolve(
         ODESYS
         , numeric
         , method = rkf45
         , maxfun = -1
         , parameters = params
         , events = [[r(t) = landing_r, halt]]
       ):

# Specialize the solution for some given value of rho, here rho=0.1

ans(parameters=[0.1]);
                  [rho = .1]

# initialize the solution procedure for some time a priori larger than
# the expected landing time:

ans(100):

# Catch the value of landing_t (meaning the value of t which fires the event)

landing_t := ans(eventfired=[1])[]
Warning, cannot evaluate the solution further right of 264.44373, event #1 triggered a halt
                   landing_t := 11.9098649021394

Next do the plots from 0 to landing_t.
For more details download the worksheets at the end of this answer.

By the way: shouldn't  g be equal to the gravitational acceleration on Mars (about 3.73 m/s² ... you seem to get a value of about 12.64 m/s²)?

Last but not least: the worksheet below provides a short analysis of the emergence of a singularity. From this analysis it is clear that the right hand side of the d𝜙(t)/dt equation takes an infinite value for some finite time and that comes from (cos(𝛾(t)) = 0 at this same time (cos(𝛾(t)) is present at the denominator of d𝜙(t)/dt equation).
This means the flight path angle 𝛾(t) is equal to -𝜋/2 at finite time (vertical drop)... which should not be possible when there is friction.
My advice: check your equations, data and units.
singularity_analysis_sand15.mw   see also   capturing_cos(gamma(t))=0.mw

I guess someone else could interpret your question differently and provide another answer.
Whatever here is mine: My_interpretation.mw (file can't be uploaded)

-0.06888275721

You write "I want to choose A to obtain the smallest absolute value of Em, preferably ~ -0.00005".
I suppose you want to say "I want to choose A to obtain the smallest absolute value of Em, preferably ~ +0.00005", don't you?

Whatever, Em is strctly negative for A > 0 and its limit when A goes to 0 from the right is -0.06888275721, which means that your requirement cannot be met.

Note that for A = 0 we have  E(𝜐) = C𝜐² + -0.06888275721 which immediatly gives the value of C... but the fitted model is unacceptable has shown in the attached file.

At the very end I'm leaning to consider your problem is not correctly posed and I suggest you to consider carefully what you want to do.
ML_sand15.mw

for i2 varying from 0 to 1

doen't nmean anything at all as the 
Error, controlling variable of for loop must be a name or sequence of 2 names
says: between for and from you must have only one name (here you have i2 and varying).

Given the past worksheets you posted I don't even understand how you can have written  (sorry if I feel that rude) such a stupid  thing?
A few correct syntaxes are

for i2 from 0 to 1 do
...
end do:

for i2 from 0 to 1 by 1/3 do
...
end do:

# etc

Note that the first syntax seems quite stupid to me in your context:

# with K2 = 0.1968052257

for i2 from 0 to 1 do
  if i2 > K2 then 
    do_something
  else
    do_something_else
  endif:
end do:

because the counter (here i2) being an integer unless "by some_value" is specified, it takes only the values i1=0 and i2=1. So a simpler operation is 

Do_this := i2 -> piecewise(i2 > K2, do_something, do_something_else)

Nothing of what you try to achieve is clear to me.
Moreover, if you are not capable to manage correctly the 2D input document mode, just as I can't neither, use the 1D input worksheet mode (your file contains a lot of breaks which generate errors (yout for anf if structures do not end correctly).

Whatever... look to the attached file to see how to write something error free A_no_error_worksheet_sand15.mw and try to use it contents to formulate correctly your problem (which I repeat I do not understand).

(done with Maple 2015 but should work the same with Maple 2018)
One_way_sand15.mw   (read carefully the last comment in thiw worksheet)

Proceed the same way for the 3D plots by replacing p1_final by p1 in C11 and p2_final by p2 in C12.

Note: If you do some browsing on this same site you should be able to find several questions about the same subject.
@acer already gave a lot of answers too [moderator: eg. here, here, here]
[moderator: If you need to export the 3D plot with colorbar in such an old Maple version then see also here.]

[moderator: colorbars received direct support for 3D plotting commands in Maple 2024.]

restart;
n := 0:
L := []:

for a from 3 to 100 do
    for b from 3 to a do
        c2 := a^2 - a*b + b^2;
        c := isqrt(c2);
        if c^2 = c2 then
            if c < a + b and a < b + c and b < c + a then
                if igcd(a, b, c) = 1 then
                    x2 := (-a^2 + b^2 + c^2)/2;
                    y2 := (a^2 - b^2 + c^2)/2;
                    z2 := (a^2 + b^2 - c^2)/2;
                    if 0 < x2 and 0 < y2 and 0 < z2 then
                        # One way
                        S := [sqrt(x2), sqrt(y2), sqrt(z2)];
                        S := S[sort(evalf(S), output=permutation)];
                        x := S[1]; 
                        y := S[2]; 
                        z := S[3];
                        n := n + 1;
                        L := [op(L), [x, y, z, sqrt(x^2 + y^2), sqrt(y^2 + z^2), sqrt(x^2 + z^2)]];
                    end if;
                end if;
            end if;
        end if;
    end do;
end do;

n;
L;

Why was sort([sqrt(x2), sqrt(y2), sqrt(z2)]) ineffective?
Look at this

S0 := [6*sqrt(110), 2*sqrt(610), 3*sqrt(649)];
sort(S0, `>`)
Error, invalid input: a numeric list is expected for sort method ``>``

# To understand this "failure" look at this 
if S0[1] > S0[2] then "yes" else "no" end if;
Error, cannot determine if this expression is true or false: 2*610^(1/2) < 6*110^(1/2)

So sort, which vasically uses an if structure, simply doesn't know how to compare these two non rational numbers.

Note that is can do it:

is(S0[1] > S0[2]);
                             true
is(S0[1] < S0[2]);
                             false

(strangely it takes a long time to be processed)

Download simplify_radical_sand15.mw

From the code below you should be able to define any binary operator you want.

Download Special_Operators.mw

AA := value( eval(A, {Omega = (t -> 1), chi = (t -> 1), g = (t -> 1), i = 1, p = (t -> 1), i = 1}) );
has(AA, {int, Int})
                             false

h-pde-M_sand15.mw

t-plot_sand15.mw

Some functions, like 'sin' for instance, leads to an unevaluated integral.

Even using evalf/Int doesn't help for the the numeric values are Float(undefined).

Download DGL_test_sand15.mw