If you use the options timestep and spacestep to pdsolve(..., numeric) to set those values to the same values used in your custom-coded solution, then the difference of the two curves will be much reduced. So, starting with 

pds:= pdsolve(
    PDE, IBC, numeric,
    time = t, range = 0 .. 1, timestep = 0.01, spacestep = 0.1
):

the final plot is 

You can see the contents as a list by

[seq](u[1..11, 1]);

You can plot it by

plot(<<seq(1..11)> | u[1..11, 1]>);

Like this:

AllWords:= (ABC::string, LW::nonnegint)->
    (op@StringTools:-Generate)~([$1..LW], ABC)
:
W:= CodeTools:-Usage(AllWords("abc", 12)):
memory used=92.01MiB, alloc change=75.49MiB, 
cpu time=719.00ms, real time=585.00ms, gc time=328.12ms

nops(W);
                             797160

Like this:

ListCoeffs:= proc(P, V::list(name))
local C, T, t;
    C:= coeffs(P,V,T);
    [seq]([degree~(lhs(t), V), rhs(t)], t= ([T]=~ [C]))
end proc:
p:= c1*x^3*y^2 + c2*x^4*y^3;
ListCoeffs(p, [x,y]);
                  [[[4, 3], c2], [[3, 2], c1]]

I answered this Question for you several months ago. Was there any problem with that Answer?

Joe's hack works by changing something before the plot is created. My hack works by altering the plot after it's created. This is used after the DrawGraph:

F:= [indets(GraphTheory:-WeightMatrix(G), fraction)[]];
P2:= subs(sprintf~("%0.3g", F)=~ sprintf~("%a", F), P); 
plottools:-rotate(P2, Pi/4);

Define it like this: 

f:= (x,y)-> 2*x^2 + 5*y^3 + 5;

The syntax f(...):= ... is unreliable, even for a function of one variable.

Your parentheses are unbalanced. Your first left parenthesis on the 4th line after unapply should be removed.

The syntax required to use color to represent the values of a function defined over a surface is quite simple, and I show it below. However, your function w (regardless of the coordinate system by which we interpret theta) is so extreme and erratic that I find it impossible to achieve visually meaningful results. This is why it has taken me so long to answer.

I tried many monotonic and compactifying transformations to w, and they did make some small improvement. But the presence of sinh with that large coefficient on x causes an extreme range of values, -10^11..10^9; and presence of cos with that large coefficient on theta causes the middle "reasonable" values to be extremely erratic with respect to theta.

So, anyway, here's the syntax to do it, and the results would look more meaningful if w were a reasonable function:

restart
:
#The color function: x is the Cartesian coordinate;
#theta is either (1) the (spherical/cylindrical)-coordinate azimuth or
#(2) the spherical-coordinate polar angle. I show option 2 below. 
#For option (1), simply change Pi/2-arctan(z(r)/r) to phi. 
#In either case, the range of w is extreme.

w:= subs(
    [x= r*cos(phi), theta= Pi/2-arctan(z(r)/r)],
    0.01503546462*(sin-sinh)(-2.365 + 9.46*x)*cos(6*theta)
):
#The surface:
(a, R, z):= (2, -8, r-> sqrt(R^2 - (r - a + R)^2)):
domain:= phi= -Pi..Pi, r= 2..3, coords= cylindrical:
plot3d([[r, phi, z(r)], [r, phi, -z(r)]], domain);

# **If** w were not such an extreme function, the following syntax would
#be all that's needed to map w onto the surface. But the prevalence of 
#extreme and erratic values makes this unsatisfactory.

plot3d(
    [[r, phi, z(r)], [r, phi, -z(r)]], domain, 
    color= [w, subs(z= -z, w)]
);

​​​​​If you got rid of the extreme coefficients of w and instead made it

w:= (sin(x) - 0.1*sinh(x))*cos(theta)

and applied the same coordinate transformation as above, then you'd get this plot:

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);