Although VV's method is fast and exact, a further substantial time improvement can be made by integerizing the matrix of rational polynomials. This is done by multiplying it by its lowest common denominator. Also, the technique below can return the LU factorization of A, which'll be beneficial if you'd like to use different bs with the same A.

#This procedure returns a factorization of A, A = d*A1, where d 
#is rational and A1 is a matrix of polynomials with integer 
#coefficients.
IntPolyMatrix:= proc(A::rtable)
local
    x, 
    AR:= convert(A, rational),
    d:= ilcm(seq(x, x= denom~(AR))),
    V:= indets(A, name)
;
    1/d, rtable(d*AR, datatype= `if`(V={}, integer, polynom(integer, V)))
end proc
:
#This procedure returns a factorization of A, A = d*P.L.U, where d is
#rational and P, L, and U are the standard matrices returned by
#LUDecomposition.
LUIntPolyMatrix:= proc(A::Matrix)
local A_denom, A_int;
    (A_denom, A_int):= IntPolyMatrix(A);
    (    
        A_denom, 
        LinearAlgebra:-LUDecomposition(
            A_int, 
            ':-method'= ':-FractionFree', 
            ':-output'= [':-P', ':-L', ':-U']
        )
    )
end proc
:
#This procedure returns the exact matrix/vector X such that A.X = b.
#X is a matrix of rational functions with integer coefficients.
SolveIntPolyMatrix:= proc(A::Matrix, b::{Vector, Matrix})
local b_denom, b_int, dPLU;
    dPLU:= LUIntPolyMatrix(A);
    (b_denom, b_int):= IntPolyMatrix(b); 
    b_denom/dPLU[1]*LinearAlgebra:-LinearSolve([dPLU[2..]], b_int)
end proc
:
#Time test (using the A and b from your worksheet):
X:= CodeTools:-Usage(SolveIntPolyMatrix(A,b)):
memory used=277.44MiB, alloc change=0 bytes, 
cpu time=1.47s, real time=1.36s, gc time=250.00ms

#Verify correctness:
[seq](x, x= simplify(convert(A, rational).X - convert(b, rational)));
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 0, 0, 0, 0]

Yes, it's easy. To do it completely, there are 2-1/2 or 3 steps:

1. Use the command LibraryTools:-Save(x, y, z, ..., filename) to create a library named filename and add x, y, and z to it. The same command can be used to add new names to an existing library. Anything that can be assigned to a name can be saved; they don't need to be procedures.

2. The global variable libname determines which libraries are visible and the order in which they are accessed (very important if the same name is used in multiple libraries). libname contains the names of the directories that Maple will look in for the libraries.

3. To make it fully automatic as you describe, you must update the value of libname in your initialization file. Make sure that you update it, not just set it. For example,

libname:= libname, "the new directory":
or
libname:= "the new directory", libname:

depending on the access order you want in case of duplicate names.

I will assume that this is an initial value problem (IVP). If it's a boundary value problem (BVP), then the code below will need to be modified a bit. It's definitely an IVP or a BVP, but which one it is cannot be determined from the text snippet provided.

If you can supply numeric values for the 4 initial conditions and 2 parameters, then the rest is easy, and done below. I just made up some numbers for those 6 values. The key commands are
dsol:= dsolve(..., numeric, parameters= [...])
and
plots:-odeplot(dsol, ...)

restart
:
ODEs:= <
    #1.85:
    diff(chi(t),t$2) + 3*H*diff(chi(t),t)       
        + (exp(-sqrt(8/3)*chi(t))*phi(t) 
            - exp(-sqrt(2/3)*chi(t))*(phi(t)-zeta*diff(phi(t),t$2)^2)
          )/sqrt(24),

    #1.86:
    zeta*(diff(phi(t),t$2)+diff(phi(t),t)/3*(9*H-sqrt(6)*diff(chi(t),t)))
        - exp(-sqrt(2/3)*chi(t))/4 + 1/2
>; 
#Parameterized initial conditions:
ICs:= phi(0)=phi0, D(phi)(0)=dphi0, chi(0)=chi0, D(chi)(0)=dchi0
:
#Extra expressions:
Extra:= <
    #1.87:
    Omega__m = 
        1 -
        (2*diff(chi(t),t)^2 + 
            exp(-sqrt(2/3)*chi(t)) *
                (diff(phi(t),t)^2 + phi(t)*(1-exp(-sqrt(2/3)*chi(t))/2))
        )/12/H^2,

    #1.88:
    q = 
        -1 - 3*Omega__m/2 - diff(chi(t),t$2)^2/2
            - zeta/4*exp(-sqrt(2/3)*chi(t))*diff(phi(t),t)^2,

    #1.89:
    omega__DE = 
        (diff(chi(t),t)^2 - exp(-sqrt(2/3)*chi(t))*diff(phi(t),t)^2)
            /3/Omega__m - 1,

    #1.90:
    Omega__DE = 
        (2*diff(chi(t),t)^2 + 
            exp(-sqrt(2/3)*chi(t)) *
                (diff(phi(t),t)^2 +
                    phi(t)/H^2*(1-exp(-sqrt(2/3)*chi(t))/2)
                )
        )/12
>;      
Params:= [H, zeta, phi0, dphi0, chi0, dchi0]:
dsol:= dsolve({seq(ODEs), ICs}, numeric, parameters= Params):
#I just made up some values:
Param_vals:= Params=~ [.1, .2, .3, .4, .5, .6];
 Param_vals := [H = 0.1, zeta = 0.2, phi0 = 0.3, dphi0 = 0.4, 
   chi0 = 0.5, dchi0 = 0.6]

dsol(parameters= Param_vals):
plots:-odeplot(
    dsol, 
    eval[recurse](
        `[]`~(t, [seq](rhs~(Extra))),
        {Param_vals[], seq(Extra)}
    ), 
    t= 0..2,
    legend= [seq](lhs~(Extra))
);

It looks like you're trying to compute the stationary vector s of the Markov chain whose transition matrix is A. For a Markov chain, the matrix Id-A will always be singular. As you've noted, you need to add a side relation; in this case, the appropriate side relation is that the sum of the entries of s is 1. Once that is done, you could use a Moore-Penrose inverse to get s, but that's not the best way. 

Below, I show three ways: 1) Least squares solution of linear system, 2) Moore-Penrose, 3) Eigenvector. In this particular case of a very small A, they produce identical results (i.e., identical at the default precision Digits=10). But for larger A, I think that 2) and 3) will take more computation time and they'll have more round-off error.

In the code below, I use ^+ as the transpose operator. Depending on your Maple version and your input mode (1D or 2D), you may need to change that to ^%T.

restart;
LA:= LinearAlgebra:
# Convenient plaintext vector display format:
Ans:= (v::uneval)-> 
    printf("\t\t%a = <%9.7g>\n", v, simplify(eval(v), zero)
):

#Your transition matrix:
A:= <0.5661180126, 0.4338819876; 0.8316071431, 0.1683928571>;
n:= upperbound(A)[1]; #matrix size
#
# Method 1: Least squares solution of linear system
#
Id:= rtable(identity, (1..n)$2, subtype= Matrix);
A1:= <Id-A | <(1$n)>>^+; #include coeffs of side relation.
R:= <(0$n), 1>; #right-side vector, including side relation.
s:= LA:-LeastSquares(A1, R)^+:  Ans(s); 
             s = <0.6571429 0.3428571>

#Check:
Ans(s.A); Ans(s.A-s);
             s . A = <0.6571429 0.3428571>
             s . A-s = <1.000001e-10 1.000000e-10>

#
# Method 2: Moore-Penrose
#

#You don't need "method= pseudo" if the matrix isn't square.
s1:= (LA:-MatrixInverse(A1).R)^+:  Ans(s1);
             s1 = <0.6571429 0.3428571>

#Compare:
Ans(s1-s);
             s1-s = <1.110223e-16 -1.110223e-16>

#
# Method 3: 
#     The stationary vector is the dominant left eigenvector of the 
#     transition matrix normalized in the 1-norm
#
s2:= LA:-Normalize(LA:-Eigenvectors(A^+)[2][..,1], 1)^+:  Ans(s2);
             s2 = <0.6571429 0.3428571>

#Compare:
Ans(s2-s);
             s2-s = <-8.689216e-11 -7.542467e-12>

#That s2 contains complex values with 0 imaginary parts that could be 
#easily removed if that's desired.

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.