Once your expression F is correctly entered, the following command is all that you need to compute s:

fsolve(F, s= 0..1);
                          0.9976968678

You did correct all the instances of ln in F. But there were many other multiplication signs missing; I added so many that I lost count. I'm not sure that I entered all of them correctly, so you should check my edited F. I also removed a huge number of unnecessary parentheses.

F:= (1-x)*R*(1/2*((1-rho*(1-s))*ln(1-rho*(1-s))+rho*(1-s)*
ln(rho*(1-s))+(1-rho*(1+s))*ln(1-rho*(1+s)) + rho*(1+s)*
ln(rho*(1+s))))+(1-x)*R*T/2*((1-s)*ln((1-rho*(1-s))*rho*
(1-s))+(1+s)*ln(rho*(1+s)*(1-rho*(1+s)))+(1-x)*(1-s^2)*
((1-4*rho+3*rho^2)*mu+(2*rho-3*rho^2)*nu))*
0.1801539058*T^(3/2)+((1-x)*R*T*rho*ln(((1+s)*
(1-rho*(1-s)))/((1-s)*(1-rho*(1+s))))-
(2*s*rho*(1-rho)*(1-x)*((1-rho)*mu+rho*nu)))* 
((-(1-x)*R*T/2*((1-s)*ln((1-rho*(1-s))*rho(1-s))+
(1+s)*ln((rho*(1+s)*(1-rho*(1+s)))+
(1-x)*(1-s^2)*((1-4*rho+3*rho^2)*mu+(2*rho-3*rho^2)*nu))*
0.1801539058*T^(3/2)-(1-x)*R*(1/2*((1-rho*(1-s))*
ln(1-rho*(1-s))+rho*(1-s)*ln(rho*(1-s))+(1-rho*(1+s))*
ln(1-rho*(1+s))+rho*(1+s)*ln(rho*(1+s))))))/
(((1-x)*R*T/2*rho*ln(((1+s)*(1-rho*(1-s)))/
((1-s)*(1-rho*(1+s))))-(2*s*rho*(1-rho)*(1-x)*((1-rho)*mu+
rho*nu)))));

If your goal for saving the data is just to re-use it in Maple, then it's a trivial one-liner to save anything created in Maple to a file. And it's an even more trivial command to import that data into a fresh Maple session. Thus, you'd might as well save the entire dsolve solution. Both the saving (command save) and importing (command read) are shown in the code below.

I also want to show you easier ways to do your entire project. Specifically, the entire thing can be done with a single call to dsolve by using the parameters option:

restart
:
H:= {h||(1..3)(z)};
H2:= combinat:-choose(H,2);
eq1:= add(mul~(H2)) - 3*(1+z)^3*(_Omega_d + _Omega_r*(1+z)) = 0; 
eq||(2..4):= seq(
    add(f)^2 - mul(f) - (1+z)*(add(f*~diff(f,z))/3 - _Omega_r*(1+z)^3) = 0,
    f= H2
);

n:= 10:
h00:= Array(0..n, i-> i/60, datatype= hfloat):
h10:= 1.-~h00: 
h20:= 1.-~.7*h00:
h30:= 3.-~h10-~h20:
ans:= dsolve(
    {eq||(2..4), seq(h||i(0)= _h||i, i= 1..3)}, numeric,
    parameters= [_Omega_r, _h||(1..3)]
):
Font:= [Times, roman, 20]:
Opts:= 
    font= Font, axes= boxed, 
    legendstyle= [location= right, font= Font], size= [1500, 1000],
    labelfont= Font, thickness= 3
:
Omega_r:= 1e-5:
h1plot:= plots:-display(
    [seq](
        (
            ans(parameters= [Omega_r, h10[i], h20[i], h30[i]]),
            plots:-odeplot(
                ans, [z, h1(z)], z= -0.9..10, 
                legend= typeset(h__1(0)= nprintf("%5.3f", h10[i])),
                color= COLOR(HUE, .85*i/n), linestyle= 1+irem(i,7)
            )
        )[2],
        i= 0..n
    ),
    Opts, 
    title= "z versus h1(z) for different initial conditions", 
    labels= ["z", "h1(z)"]
);

In Maple 2019 or later, with 1D input, this syntax can be used to get the same thing:

h1plot:= plots:-display(
    [
        for i from 0 to n do
            ans(parameters= [Omega_r, h10[i], h20[i], h30[i]]);
            plots:-odeplot(
                ans, [z, h1(z)], z= -0.9..10, 
                legend= typeset(h__1(0)= nprintf("%5.3f", h10[i])),
                color= COLOR(HUE, .85*i/n), linestyle= 1+irem(i,7)
            )
        od
    ],
    Opts, 
    title= "z versus h1(z) for different initial conditions", 
    labels= ["z", "h1(z)"]
);

Now here's how to make your second plot:

i:= 2:
ans(parameters= [Omega_r, h10[2], h20[2], h30[2]]):
allplot1:= plots:-odeplot(
    ans, `[]`~(z, [h||(1..3)(z)]), z= -0.8..8,
    color= [red, blue, green],
    legend= [seq]( 
        typeset(h__||j(z), " with ", h__||j(0)= nprintf("%5.3f", h||j||0[i])), 
        j= 1..3
    ),
    Opts,
    title= typeset(h(z), " vs. ", z, " (case ", i, ")"), labels= [z, h(z)]
);

To save everything, all you need to do is

#The filename must end with .m!
save Omega_r, h10, h20, h30, eq||(2..4), Opts, ans, "zvshGR.m":

Now I start a new Maple session and create the plot that you mentioned in Question but didn't create in the worksheet:

restart:
read "zvshGR.m":
i:= 5:
ans(parameters= [Omega_r, h10[i], h20[i], h30[i]]):
#We need to *algebraically* solve the original ODEs for their highest-order 
#derivatives:
DS:= eval(solve({eq||(2..4)}, diff({h||(1..3)}(z), z)), _Omega_r= Omega_r): 
plots:-odeplot(
    ans, [z, eval(((add@rcurry(diff, z))/add)([h||(1..3)(z)]), DS)], z= -0.8..8
); 

A .lib file, which must be used with a matching .ind file, is a very old forerunner of the current .mla file. I think that the best thing to do would be to use the command Library:-ConvertVersion to convert it to a .mla file and then upload that to the Cloud. You may need to do the conversion with a local copy of the .lib.

You can increase the order as needed. This can be done automatically by this little procedure:

MySeries:= proc(F::algebraic, X::{name, name= algebraic}, order:= Order)
local S:= series(F,X,1+order), e:= op(2,S);
    `if`(e < 0, series(F, X, 1+order-e), S)  
end proc: 

Try this:

Grid:-Set(fun, UtilsIORelation, 'model', 'vars', 'conds');

If that doesn't work, try

Grid:-Set(fun, UtilsIORelation:-findMatchingStructures, 'model', 'vars', 'conds');

Like this:

eval([x,y], solve({y=y2}));

Like this:

select(isprime, 2^~[$1..100] -~ 1);
      [3, 7, 31, 127, 8191, 131071, 524287, 2147483647,  2305843009213693951, 618970019642690137449562111]

There are several ways to improve the efficiency of this, and primes of the form 2^n-1 are widely studied (they're called Mersenne primes). The most basic improvement is that the exponent must itself be prime because 2^(a*b)-1 is divisible by 2^a-1 and 2^b-1. Thus, the sequence above can be generated by

select(isprime, 2^~select(isprime, [$1..100]) -~ 1);

You can do, initially,

phi:= `#mo("phi")`;

and then use phi as normal. You could do this with any symbol.

I don't think that this has anything to do with palettes; although, I never use palettes, so I'm unsure about that. I think that any symbol that you see on a palette, as well as nearly any symbol or combination of symbols that you can think of (mathematical or otherwise), can be created via a mechanism akin to the one shown above.

I don't think that it makes sense to integrate directly with respect to a unit, but you can include units in limits of integration. So, you can do this:

with(Units:-Standard):
Intat(10*Unit(m/s^2), t= Unit(s));
combine(value(%), units);

or this:

Int(10*Unit(m/s^2), t= 0..Unit(s));
combine(value(%), units);

Intat is something between a definite and an indefinite integral (antiderivative): It evaluates the antiderivative at a single value of the integration variable.
 

Edit: I just read your followup to Rouben that you posted while I was writing this, so I understand now that you want to do more than just integrate this function. Getting your worksheet now.

Original answer: The function is piecewise polynomial with exact rational coefficients and breakpoints. It's trivial to integrate it symbolically (even if it weren't smooth or continuous) and get the exact value.

int(Z(x), x= 5/3..5);

             361954048721879 / 361954048856400

Approximating the function by some other function is a bad idea. The computational effort required to get an approximating function is far greater than that required for integration.

You have defined

G:= array(1..N, 1..N);

You should change that to

G:= Array(1..N, 1..N);

Then G will print. The command array is deprecated; you should never use it.

Here is a much more general way by which we can work with differential operators whose behavior mimics that of D itself, and when applied to arguments, the expressions are converted to diff form. I'n going to call my operator DiffOp, but feel free to change that name.

`evalapply/DiffOp`:= proc(DO, Y)
    local y,x;
    if procname::indexed then
        y:= op(procname); x:= Y[];
        evalindets(
            evalindets(
                op(DO),
                {
                    'D'(identical(y)), 
                    And(
                        typefunc(
                            identical(y), '`@@`'(identical(D), anything)
                        ),
                        function(identical(y))
                    )
                },
                d-> convert(d(x), diff) 
            ),
            identical(y),
            f-> f(x)
        )            
    else
        y:= Y[];
        DiffOp[y](
            subs(
                {D@@0, `1`}=~ y,
                evalindets(
                    evalindets(
                        subs(D= D(y), op(DO)),
                        '`@@`'('D'(identical(y)), anything),
                        d-> (D@@op(2,d))(y)
                    ),
                    {procedure, And(appliable, `module`)},
                    p-> p(y)
                )
            )
        )
    fi
end proc
:

#Example: 
# Use D@@n instead of D^n. The @@ is Maple's operator for functional iteration. 
# D^n will be interpretted as the nth algebraic power of the first derivative.
# Use `1` or D@@0 (the choice is yours) for the identity operator (which will become
# y(x) when the differential operator is applied.

did:= x*D@@3 + 144*sin*D@@12 + 7*D + `1`;
This output makes more sense when prettyprinted in a worksheet:
        did := x @@(D, 3) + 144 sin @@(D, 12) + 7 D + 1

DiffOp(did)(y)(x);

  /  3      \                 /  12      \                      
  | d       |                 | d        |     / d      \       
x |---- y(x)| + 144 sin(y(x)) |----- y(x)| + 7 |--- y(x)| + y(x)
  |   3     |                 |   12     |     \ dx     /       
  \ dx      /                 \ dx       /                      

                 

I will explain to you why this doesn't work with eval or subs, and I'll even teach you an alternative that does work, but I must emphasize that your whole idea along these lines is very, very, VERY cruddy and prone to error; so, you must never do this or ask a Question about it again. Your subsequent Question, already Answered by Kitonum, shows that you're well on your way to having an alternative to this mess.

Why it doesn't work: The issue is order of evaluation. An expression being modified by subs or eval is evaluated before the substitutions are made. (This is simply the default order of evaluation for the vast majority of Maple commands.) So, at the time of evaluation, i has no value. That doesn't make sense to op. There are many Maple commands that are prepared to handle situations where a parameter that is usually numeric is passed symbolic (i..e., has no value). And there are many other commands, such as op, that are not prepared to handle that.

An alternative that works: Use a procedure to change the order of evaluation:

MyC:= i-> coeff(op(i,did), D, op(2,op(2,op(i,did))));
MyC(3);

Even better:

MyC:= (i::nonnegint)-> coeff(op(i,did), D, op([i,2,2], did));
MyC(3);

Note that nested op calls can be compressed by using a list as the first argument.

Use an embedded if statement rather than `if`, like this:

str:= "A";
x:= 10;
str:= sprintf(
    "%s it was%s 10", str, if x++=10 then "" else x:= 9; " not" fi
); 

The benefit of an embedded if statement is that any number of statements can appear in its then and else clauses. It's only the value of the final of these statements that becomes the return value of the embedded statement.

Your computation can be done by very simple arithmetic in compiled code. For the n=8 case, that avoids (2^8)^3 ~ 16 million calls to Bits:-Split, (2^8)^2 + (2^8)^3 calls to ListTools:-Occurences, and the creation and garbage collection of the same number of lists. My code below does the n=8 case in 17 seconds total on one processor. Further gains can be obtained by

• Writing the code directly in C so that you can use whole-word bit operators.
• This code can be parallelized with Grid or Threads (Threads is usually faster),
• Perhaps some speedup can be achieved by using 2-byte integers directly in C.

I made a great number of changes with the arithmetic. Please verify the accuracy of this by comparing against known-correct code using the same random vector M.

restart
:
BuildA_M:= proc(
    A::Array(datatype= integer[4]), M::Array(datatype= integer[4]),
    r::integer[4]
)::integer[4];
local 
    i0::integer[4], i::integer[4], j0::integer[4], j::integer[4],
    m::integer[4], t0::integer[4], t::integer[4],
    p::integer[4], S::integer[4]
;
    for i0 from 0 to r-1 do
        for j0 from 0 to r-1 do
            S:= 0;
            for t0 from 0 to r-1 do
                p:= 1; i:= i0; j:= j0; t:= t0; m:= M[t0+1];
                while i*t <> -j*m do            
                    if (i*t - j*m) mod 2 <> 0 then p:= 1-p fi;
                    i:= i/2; t:= t/2; m:= m/2; j:= j/2
                od;
                S:= S+p
            od;
            A[i0+1, j0+1]:= S
        od
    od;
    0
end proc
:
BuildA:= Compiler:-Compile(BuildA_M):

LArip:= proc(n::posint)
local 
    r:= 2^n, 
    M:= Array(1..r, combinat:-randperm(r) -~ 1, datatype= integer[4]),
    A:= Array((1..r)$2, 0, datatype= integer[4])
;
    BuildA(A, M, r);
    (A,M)
end proc 
:
CodeTools:-Usage(LArip(8));
memory used=263.30KiB, alloc change=0 bytes, 
cpu time=16.55s, real time=16.68s, gc time=0ns

​​​​​​