@Teep 

From Optimization/General/Options help page:

depthlimit = posint
Set the maximum depth of the branch-and-bound tree. This option is only used by the Optimization[LPSolve] command for integer programs. 
The default is the maximum of 2 and ceil(3*n/2).
The internal workspace allocated by the integer programming solver depends on the value of depthlimit.  If, during the computation, the workspace is found to be insufficient, an error is issued indicating that depthlimit must be increased.

The attached worksheet presents some tracks on a simplified problem (5 jobs and 5 machines):

1. For this n=26 variable problem the default value of depthlimit is 114 and wasn't large enougth to avoid the

   Error, (in Optimization:-LPSolve) maximal depth, 500, of branch and bound search is too small; use depthlimit option to increase depth

2. Using depthlimit=1000 returns a solution with performances

   memory used=3.47MiB, alloc change=1.93MiB, cpu time=3.70m, real time=3.71m, gc time=0ns

3. Using depthlimit=600 returns the same solution with better performances

   memory used=2.21MiB, alloc change=1.62MiB, cpu time=5.10m, real time=5.13m, gc time=0ns
   Which suggests to use only the necessary depthlimit value, no more (which is rethorical)

4. Using depthlimit=500 returns the same error as with the default depthlimit value

5. Using again the same depthlimit=500 value but restricting the number of nodes to explore (nodelimit=10^5)
   gives almost the same solution than depthlimit=600 with nodelimit=0 ... but it's  about seven times faster:
   memory used=1.97MiB, alloc change=1.55MiB, cpu time=43.61s, real time=43.67s, gc time=0ns

So I suggest you not to keep increase depthlimit (which penalizes the computational time), but set it to  a value reasonably larger than the default one (182) and increase progressively nodelimit until the  optimal solution no longer evolves.

5 jobs and 5 machines  Some_tracks_5_5.mw   !!! read red text below
5 jobs and 8 machines  Some_tracks_5_8.mw   !!! read red text below

When opening again those same files with Maple 2015 (the version I use) suprious errors appear due to some invisible character in some command blocks.
Here is the (partially) rewritten code for the 5 jobs and 8 machines case 
Some_tracks_5_8_(new).mw

XT_sand15.mw provides a very simple solution.

Nevertheless changing the x range or the tvalue will very likely make this error to appear

Error, (in dsolve/numeric/bvp) Newton iteration is not converging

Despite it is a very classical error for which some recipes do exist to fix it, see dsolve/numeric_bvp/advanced help page and the so-called continuation method, it can be hard to find the correct way to use this trick.
I blindly tried some possibilities but even I can't say any of them gave satisfaction.
Another possibility might be ti initoiate the dsolve process by giving some approximate solution  (search for approsoln in dsolve/numeric/BVP help page).

Generaly this kind of issue cannot be fixed by applying blindly those recpes and the one which is best placed to use them in an enlightened manner is the one who understands the physics underlying the differential system to be solved, i.e., you.

When several quantities have the same range of variation and this range is much greater than the relative variations between these same quantities, there is no choice but to represent the absolute or relative differences in order to highlight the differences.

Q_SEPARATE_sand15.mw

(if this suits you change the vertical label accordingly)

You wrote "i have four equations i want make them be bigger or smaller than zero" which does mean anything.
Do you mean "I want their solutions to be larger or smaller than 0"?
More precisely doyou mean "I want all w, p, chi and beta[1] to be positive or all w, p, chi and beta[1] to be negative" ?

If it is so there is no 6-uple (S1, S2, S3, T1, beta[2], k) such that all w, p, chi and beta[1] have the same sign.
It is indeed very simple toprove that chi*p is always negative, meaning chi and p have different signs:
 

Download condition-parameters_sand15.mw


If you want to find a 6-uple (chi, w, p, k, beta[1], beta[2]) such that botg S1, S2, S3, T1 are positive look here condition-parameters_2_sand15.mw
 

First point: avoid using multiple restart command in your worksheet.

In my answer I deleted all which is before the second restart and focused on this error

p1 := contourplot(P, u=-1..1, v=-1..1, scaling=constrained,
    colorscheme="DivergeRainbow", contours=[seq](x, x=-0.5..0.5,0.05)):
Error, (in plot/iplot2d:-Levels) could not evaluate expression

When you get such a message the first thing to do is to display tha quantity to plot (P) to try and see if it is plottable.

This is not the case in your worksheet because P depends on 6 parameters k1, s1, ..., s5 beyond u and v.
So give them all numeric values and and plot. Illustration below:

Download Why_cannot_P_be_plotted.mw

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