Your general solution isn't general enough. If f is *any* differentiable function such that f(0)=1, then x0 = y0 = f is a solution to your equation. 

A basic fact about systems of equations is that the number of things that can be solved for is at most the number of equations. 

If L is the list of expressions, then all that you need to do is

R:= Grid:-Map(int, L, x= 0..1);

For numeric integration, do 

R:= Grid:-Map(int, L, x= 0..1, numeric);

The command will automatically apply some rudimentary load balancing.

It's possible to get the exact solutions for all 6 functions by direct use (i.e., without using dsolve) of the Method of Undetermined Coefficients. If you try to do this with dsolve, it will get the first 5, but the 6th took longer than I cared to wait. I killed it after more than an hour. The method below gets it in seconds.
 

Download UndeterminedCoeffs.mw

The plot command is plotting over a real interval. It can't be expected to find on its own every turn of a function, especially considering that there could be infinitely many of them (such as for the Weierstrass function^[*1]). If there are values of n that you want to make sure are included, you can use the sample option. In this case, we might as well make the sample all the integers in the interval:

f:= n-> ceil(sqrt(4*floor(n))) - floor(sqrt(2*floor(n))) - 1:

Edit: My intention was to use the OP's original function, but I mistakenly copy and pasted Tom Leslie's modification of the original. So, make that

f:= n-> ceil(sqrt(4*n)) - floor(sqrt(2*n)) - 1:
plot(f, 10...100, sample= [$10..100]);

[*1] Weierstrass function: f:= x-> Sum(cos(3^n*x)/2^n, n= 0..infinity). It is continuous everywhere and differentiable nowhere.
 

You should use D for derivatives in initial and boundary conditions. I don't know why that is, but it was the only change that I needed to make to your code. Sometimes the diff form works, but the D form always works and is a lot less to type.

sys:= {
    diff(phi(eta), eta$2) + 5.261282735*f(eta)*diff(phi(eta), eta) - 
    2.630641368*phi(eta) = 0,
 
    1.059704409*diff(theta(eta), eta$2) + 6.176017503*f(eta)*diff(theta(eta), eta) 
    + 21.03607964*diff(f(eta), eta$2) + 0.5*phi(eta) = 0,
 
    diff(f(eta), eta$4) - 1.052256547*diff(f(eta), eta)*diff(f(eta), eta$2) + 
    1.052256547*f(eta)*diff(f(eta), eta$3) + 5.165076420*diff(theta(eta), eta) + 
    5.261282735*diff(phi(eta), eta) = 0,

    D(phi)(0) = 1 + 0.5*(D@@2)(f)(0),
    D(phi)(1) = 0.5*(D@@2)(f)(1), 
    f(0) = -0.5, f(1) = 0.5, 
    phi(0) = 1, phi(1) = 0, 
    theta(0) = 1, theta(1) = 0
}:
dsol:= dsolve(sys, numeric):
plots:-odeplot(
    dsol, `[]`~(eta, [f, phi, theta](eta)), eta= 0..1, 
    legend= [f, phi, theta]
);

This is an attempt to Answer your followup Question.

A single summation with two indices, say, Sum(f(i, j), i+j= 0..N), is equivalent to the nested summation

Sum(Sum(f(i, j-i), i= 0..j), j= 0..N)

So, if you define your function as 

P:= (x,y)-> Sum(Sum(alpha[i, j-i]*x^i*y^(j-i), i= 0..j), j= 0..N)

then you can differentiate it like any usual function of x and y.

You began your worksheet by assigning numeric values to "symbols". Then you used those symbols in some formulas, and then you obtained or attempted to obtain some purely numeric solutions, such as plots. This is almost always a bad order to do things---whether you're using a symbolic language like Maple, a numeric language like Matlab, or just doing it with pen and paper: It may seem convenient at first, but it becomes inconvenient when you want to change something. And, if you're using a symbolic language, that order may prevent you from using some of the symbolic features. For example, in the worksheet below, I use some of the basic commands for manipulating symbolic polynomials. These don't work as well when the polynomials contain numeric constants with decimal points (called floats). (They still work, but not as well.)

The better approach is to start with a fully symbolic presentation of the problem and then incorporate the numerics at the end.

The central formula (equation) of your worksheet (the one in the solve command) can be easily converted to a polynomial in x. And, if you keep it symbolic, most above-average high-school algebra 2 students could do that conversion with pen and paper. If decimal numbers are used before this step, the result becomes a difficult-to-read mess. There are many benefits to using a polynomial. For example, fsolve automatically returns multiple solutions for polynomials, but not for other equations.

The worksheet below runs many, many times faster than the original approach of solving for 100 values of t.
 

Download ProcSubs2Levels.mw

If you'd like the horizontal axis to be scaled by V__line (as VV did), I can make a small adjustment to do that.

You could use the top-level longstanding galois command for the information that you want.

galois(x^5 + 20*x + 32);
              "5T2", {"5:2", "D(5)"}, "+", 10, {"(1 2 3 4 5)", "(1 4)(2 3)"}

The help page ?galois describes what all that means.

On the other hand, any place where a space is allowed, you can use multiple spaces, or line breaks, or any combination of those.

The command convert(..., binary) is not very useful for doing further arithmetic; it's mainly to display as output a binary representation of a number. An extremely efficient alternative is Bits:-Split, which returns a list of 0s and 1s. It's trivial to add a list with add (after Maple 2018, I believe) or `+`@op (in any Maple). To create your histogram, you could do

Statistics:-Histogram((`+`@op@Bits:-Split)~([$0..2^16]), binwidth= 1);

Another option that returns the same list is convert(..., base, 2). This is much slower than Bits:-Split (which is compiled code), but it's a command worth knowing because the 2 can be replaced by any positive integer other than 1; whereas Bits:-Split can only work with bases that are powers of 2. And it's not a very slow command; it's just relatively slow compared to the super-efficient Bits:-Split. The above command takes 2 seconds on my computer. If I use convert(..., base, 2), then it takes 7 seconds.

Both of these commands return the digits in order from least significant on the left to most significant on the right. Although that looks backward to us due to the way we usually write numbers, it is the form most convenient for further processing. It is only due to cultural rather than mathematical reasons that it looks weird. 

The dsolve has returned an equation whose left side is y(x) and whose right side is some expression containing x but not y. We call that an explicit solution. You don't want to carry around that left side for the remaining computations. Change the line

C := eval(N__1, x = 200)

to

C:= eval[recurse](y(x), [N__1, x= 200])

There's no need to "estimate" the distance; it can be done exactly. Modify your distance function, like this:

dist:= proc(p::[integer,integer], q::[integer,integer])
local H;
    min((add@abs~)~([seq]([p[2]-H, p[1]-q[1], H-q[2]], H= [1,10])))
end proc  
:
Edgs:= {seq}([({op}@lhs~)(e), dist(rhs~(e)[])], e= combinat:-choose(pts2, 2)):
G1:= GraphTheory:-Graph(Edgs):
(distance, Tour):= GraphTheory:-TravelingSalesman(G1);
            distance, Tour := 38, [A, D, C, B, E, A]

I'm working on a method to highlight the trail with arrows.

Your troubles were caused by a hard-to-see typo. The subexpression beta(sin(alpha)/sin(beta)) occurs in the denominator of eq1. That beta is not multiplied by its following expression. To get implied multiplication in 2D Input, you need a space after beta. Or, you can use an explicit multiplication symbol *. If I make this correction and specify variable ranges to fsolve, then I get a solution close to (but not identical to) the one you mentioned earlier.

Download CriticalPressureFsolve.mw

If you

1. use 15 Digits of precision,
2. use fsolve instead of solve,
3. use Pi instead of 3.1415, and
4. search for complex solutions,

then the imaginary parts of those solutions will be very small, and you can ignore them. After you've entered the three equations, do this:

Digits:= 15:
PI:= Pi:
solC:= fsolve((Eq:= {eq||(1..3)}), complex);
         /                        8                      -11    
solC := { P = -1.44914053377011 10  + 1.81361988528669 10    I, 
         \                                                      
                                                 -21    
  alpha = -3.21040506122157 - 2.69567382732802 10    I, 

                                                -21  \ 
  beta = -2.50058386568901 - 1.52201908062012 10    I }
                                                     / 
solR:= eval(solC, I= 0); #Remove imaginary parts
         /                        8                             
solR := { P = -1.44914053377011 10 , alpha = -3.21040506122157, 
         \                                                      

                          \ 
  beta = -2.50058386568901 }
                          / 
evalf(eval(Eq, solR)); #Check residuals
       /       -12              -12             -10     \ 
      { -6.9 10    = 0., -3.3 10    = 0., 1.1 10    = 0. }
       \                                                / 

Red text is your input; blue text is Maple's output..

The command mtaylor does it. The code below assumes that y__i is some analytic expression in the variables x[1], x[2], x[3].

V:= [x[1],x[2],x[3]]:
ord:= 4:
sort(mtaylor(y__i, V, ord+1), V, ascending);