The  modulus of solnum is the simplest case to deal with.
Given it's theoritical range of variation [0, +oo) the effective range the grid size "captures" can go from 0 to 10⁴ with a [300, 300] grid.

So it is impossible to get a smart representation of the level curves unless you lessen this range of variation. One simple way to do this is to apply the log function (or some surd function) to the expression to plot.
Using a "compressor"  such as log, for instance, introduces extra issues because the quantity to plot must be strictly positive (whichis why I concentrated on |solnum| ).
Note that a "compressor" is a trick which enables to get a smart plot, but the level curves or the colored region do not give a faithful representation.

The main idea to have a good spacing between level curves is to aggregate pieces of |solnum| contourplorts instead of |solnum| ina row.
Examples are provided in the attached file.

To reduce the size of the grids and still get a nice picture I suggest you to account for the z-periodicity of solnum (see the attached file).

At the end you will be able to get something like this by tuning the regularizer (f), the ranges of |solnum| pieces (r1, r2, r3), the contour lists or the coloring.


The plots are removed in the the attached worksheet (size 66 Mb instead)
plot-help_sand15.mw

This is not because Maple is powerfull that not using what  one does by hand is a useless strategy. Which in your case means that replacing k by a complex expression before trying to simplify this latter is a very bad way to proceed.

BTW: A syntax like  f(x) := something is not accepted in Maple 2015 (the version I use), So I used these f := x -> something notations instead. Feel free to modify them in your own code.

More of this there is no need to define Phi(x, t) the way you did because Phi(x, t) = C*Omega(x, t) where C is some constant

Download pdetest2_sand15.mw

After the M, b := ... line insert these to verify what b is:

b;


# So b[3] is not assigned:

b[3];
Error, Vector index out of range

The error comes from you writting

M, b := GenerateMatrix(eqs, [lambda[1], lambda[2]])

Which redefines b as a zero vector of dimension 2: thus  b = <0 | 0> and the b[3] term in M generates an error.

Cure: write instead

M, B := GenerateMatrix(eqs, [lambda[1], lambda[2]])

to keep b unchanged

S_sand15.mw

@nm 

As far as I can remember, exportplot, called from cmaple, has never worked properly for versions 15 through 2021, whether on Windows or Mac OS.

I always used this kind of workaround:

d := currentdir()

# f := cat(d, "/test.png »);  # doesn't work with maple 2015 (driver unknown)

f := cat(d, "/test.jpeg");

p:=plot(x^3, x = -8 .. 8, color = "blue",axis=[gridlines=[10,color="red"]]):

g := seq(plot(y, x=-8..8, color=red), y=[seq](-500..500, 100)), seq(plot([[x, -500], [x, 500]], color=red), x=[seq](-8..8, 2)):

pg := plots:-display(p, g):    

plottools:-exportplot(f, pg);  

Here is the file exported from Cmaple (Max OSX)

Note that the rendering still remains quite bad.

... which is probably non optimal and whose efficiency could likely be enhanced by someone who have enough time to spend on it.
Note that an obvious way to fasten this code is to add the numpoints option in the plot command. For instance numpoints=30 makes the code about seven times faster.

The code is generic in the sense it will enable you to

• plot any of your 6 functions
• versus any of one of the 4 parameters,
• giving to any one of the 3 remaining ones a list of values,
• and to the last 2 parameters a single value.

Generic_2d_plot_sand15.mw

The plots below are those you asked for

But as I said earlier the code enables you to draw any expression built from traceable quantities (those  listprocedure contains):


Several examples are given in the attached file, among them:

@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