If you were going to define and declare a new constant, I don't see any point in making it anything other than global. Thus there's no point in having the value of constants revert when you exit a procedure. 

But perhaps you can convince me otherwise on either of those two points. 

Okay, with your additional information, I now understand your Question, which otherwise seemed trivial.

To answer your titular question: You get solve to ignore certain independent variables by including an identity clause in the solve command. So, here is the solution to your example:

ind_V:= {r, phi, Z}: #independent variables
Eqn:= 
    8*r^5*_C1*((alpha[4,4] - alpha[5,5])*cos(2*phi) 
    + sin(2*phi)*alpha[4,5]) = 0:
Sols:= {
    solve(
        identity(Eqn, phi), 
        indets(Eqn, And(name, Not(constant))) minus ind_V
    )
};
Sols2:= remove~(evalb, Sols); #Remove trivialities.

    {{_C1 = 0}, {alpha[4,4] = alpha[5,5], alpha[4,5] = 0}}

#Optional step: Put solution in the exact and/or form that you 
#presented:
OrAnd_form:= `or`((``@`and`@op)~(Sols2)[]);

    (_C1 = 0) or (alpha[4,4] = alpha[5,5] and alpha[4,5] = 0)

An identity clause is currently limited to declaring only 1 variable as the independent variable. It is obvious in this case that that variable should be phi, not r.

[Edit: This Answer, which seems trivial, was based on the original version of the Question. The OP subsequently added much clarifying material, and a new Answer is needed.

I hereby grant permission to any Moderator to delete this Answer if they feel, as I do, that that is appropriate. I don't think that it'd be appropriate for me to delete it myself, so I recuse myself.]

You wrote: 

• I would like to solve this type of equations for all values of the variables, namely phi in the example above.

If you want to solve for phi, the command is solve(..., {phi}) (complete example below).

• If I do not choose the variables for which to solve, I get phi=phi as one of the equations in the solution.
• If i choose all the variables except phi, i get an the expression for the constants containing the variable phi.

You seem to have tried every combination except the only obvious one: Choose phi (and phi alone) as the variable for which to solve.

Eqn:= 
   sin(phi)*_C1*a[3,5] + cos(phi)*_C1*a[3,4] = 
   sin(phi)*_C11*a[1,3] - cos(phi)*_C11*a[2,3]
;
solve({Eqn}, {phi});

• Is there a way to solve these equations automatically?

I consider  the above "automatic". If you do not, we can work on it---there are many options.

• Or do I have to separate them manually using collect and coeffs?

No, you shouldn't need to use those commands just to solve equations.

• Or could the solution be the comand Parameters from Physics package?

No, and don't you think that that would be a rather obscure way to perform one of the most important functions of a computer-algebra package---solving an equation?

You're using d and dt as if they were variables when you obviously intend them to be a differential operator. Your equation can be entered like this:

dl1:= C[th]*diff(T[ovn](t), t) = 
    P[el](t) - (T[ovn](t) - T[a](t))/R[th] - A*sigma*(T[ovn](t)^4 - T[a](t)^4);

Note in particular the derivative: diff(T[ovn](t), t). The above is written in the "1D" (aka plaintext aka "Maple Input") form that NM was refering to. You can copy-and-paste these lines into a worksheet. It'll work regardless of the input mode.

I don't think that there's much hope of getting a solution as long as there are unknown functions of t (such as P[el](t) and T[a](t)) present. The dsolve will just return NULL, I think because it gives up, not because it thinks that there's no solution. (There is definitely a solution assuming that the unknown functions are suitably smooth.)

You might as well define the nonhomogenous part as 

F:= t-> (P[el](t) + T[a](t)/R[th] + A*sigma*T[a](t)^4)/C[th];

This won't help dsolve to solve it; it's just to keep things tidy.

There's no symbolic representation of the solution other than RootOf(...). But you can do

fsolve(beta = cot(beta), beta= 0..Pi/2);

Eq:= a*x+b*y+c = 0:
expand(solve({Eq}, {y})[]);
 

The following will work for any set of equations (regardless of whether 0 is on either side). Like Kitonum's solution, it puts terms with derivatives on the left side and those without derivatives on the right side.

SepDiff:= (s::set(`=`))->
    local Z, S:= (lhs-rhs)~(subs(0= Z, s));
    subs(Z= 0, map(e-> `=`((1,-1)*~selectremove(has, e, diff)), S))
: 

To use it in your case, simply do

SepDiff(sols1);

If you have a working plot3d command such as 

plo3d(u, x= -3..3, t= -3..3);

then simply change it to

plots:-contourplot(u, x= -3..3, t= -3..3);

Note that I use u rather than [u]. The [u] is acceptable to plot3d but not to contourplot.

Here is the solution with plots for the equation solved as a BVP and solved via HPM (Homotopy Perturbation Method). I didn't make the table, but you could extract the data for the table from the matrices in the plots (indeed, for 201 values of eta in 0..1) just like I did in the code below. To my mind, all that matters is the maximum error over the interval 0..1, so that's what I computed. I think that the table was just filler for a relatively short paper.
 

Download BVPvsHPM.mw

If your goal is simply to represent objects defined in real vector spaces with more than 3 dimensions, you can use animation and/or color. Animation makes time a dimension. There are several schemes (RGB, HSV, etc.) for squeezing in 3 dimensions of color. 

I can't look at your worksheet because MaplePrimes doesn't accept the .maple format. If you resave and repost it as a .mw, I'll be able to look at it. This should be no problem because your work shouldn't require any auxilliary data.

There are two errors in the printed statement of the problem:

• R and Re are used to represent the same thing. I just changed Re to R.
• The first initial condition for the third equation should be F₂(0)=0, not 
  F₀(0)=0.

Here's my solution, which matches what's stated in the paper:

Download SequentialDsolve.mw

The command is currentdir. Example:

currentdir("C:/Users/carl/desktop");

The forward slash / can always be used as the directory separator in Maple filenames, even on Windows.

A distinction needs to be made between equations and assignments. This distinction can be a slightly difficult concept to grasp for new programmers. It is often the first nontrivial concept one encounters when learning a first computer language. This distinction exists in all computer languages, but it is more pronounced in symbolic languages such as Maple, where variables can be meaningful even when they don't have values (we say that such variables are symbolic).

Equations are declarative statements (in the logical/linguistic sense). They can be true or false depending on the values of their variables. With symbolic variables, "truth" can be more subtle: An equation may be true for all, for some, or for no values of its variables. In Maple, equations are expressed with the equals sign =.

On the other hand, assignments are imperative statements (i.e., commands) that force a variable to be a certain value. (That certain value may possibly still contain symbolic variables.) Thus, there is no question of their truth. In Maple, assignments are expressed with := or occasionally with the assign or unassign commands. Assignments are often expressed in mathematical writing with "let", as in "Let x = 3."

Your problem would be best done using equations. Like Acer, I disdain making assignments to the directly meaningful variables of a problem---those variables that are meaningful outside the context of a computer program. 

Here is the solution to your problem. Note how conveniently the Solution is expressed: It shows the numeric values of all variables in one place. These numbers are not assigned to the variables. Rather, solve is saying that if the variables had those values, the equations would be true.

System:= {
   totalsalesx = 60.1 + tax + profit,
   tax = 0.3*profit,
   profit = 0.1*totalsalesx
}:
Solution:= solve(System);
           {profit = 6.908045977, tax = 2.072413793, totalsalesx = 69.08045977}

Here is a procedure for it:

DroppedContourPlot3d:= proc(
   Z::name= algebraic, #function to plot
   X::name= range(realcons), #x-axis-parameter range
   Y::name= range(realcons), #y-axis-parameter range
   {
      #fraction of empty space above and below plot (larger "below"
      #value improves view of dropped-shadow contourplot):
      zmargin::[realcons,realcons]:= [.05,0.15],
      contouropts::list({name, name= anything}):= [],
      surfaceopts::list({name, name= anything}):= []    
   }
)
option `Author: Carl Love <carl.j.love@gmail.com> 19-Nov-2019`;
local
   LX:= lhs(X), RX:= rhs(X), LY:= lhs(Y), RY:= rhs(Y),
   C:= plots:-contourplot(
      rhs(Z), X, Y,
      #These default plot options can be changed by setting 'contouropts': 
      'contours'= 5, 
      'filled', 'coloring'= ['yellow', 'orange'], 'color'= 'green', 
      contouropts[]
   ),
   P:= plot3d(
      rhs(Z), X, Y,
      #These default plot options can be changed by setting 'surfaceopts': 
      'style'= 'surfacecontour', 'contours'= 6,
      surfaceopts[]
   ),
   U, L #z-axis endpoints after margin adjustment
;
   #Stretch z-axis to include margins:
   (U,L):= ((Um,Lm,M,m)-> (M*(Lm-1)+m*Um, M*Lm+m*(Um-1)) /~ (Um+Lm-1))(
      zmargin[],
      (max,min)(op(3, indets(P, 'specfunc'('GRID'))[])) #actual z-axis range
   );
   plots:-display(  #final plot to display/return
      [
         plots:-spacecurve(
            { 
               [  #(y,z)-plane backwall:
                  [lhs(RX),rhs(RY),U], [rhs(RX),rhs(RY),U],
                  [rhs(RX),rhs(RY),L]
               ],
               [  #((x,z)-plane backwall:
                  [rhs(RX),rhs(RY),U], [rhs(RX),lhs(RY),U],
                  [rhs(RX),lhs(RY),L]
               ] 
            },
            'color'= 'grey', 'thickness'= 0
         ),
         plottools:-transform((x,y)-> [x,y,L])(C), #dropped contours
         P
      ],
      #These default plot options can be changed simply by putting the option
      #in the DroppedCountourPlot3d call:
      'view'= ['DEFAULT', 'DEFAULT', L..U], 'orientation'= [-135, 75], 
      'axes'= 'frame',
      'labels'= [lhs(X), lhs(Y), 'z'], 
      'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'axesfont'= ['TIMES', 12],
      'projection'= 2/3,   
      _rest
   )
end proc
: 

Example usage:

DroppedContourPlot3d(z=x^2+y^2, x= -2..2, y= -2..2);

I didn't plot the same function that you specified because it's too difficult to read your unformatted Question, but it should be straightforward to plot it as above.

My procedure automatically chooses a plane for the 2d plot to optimize the view; it doesn't necessarily use plane z=0.

Since you tried partial fractions, I thought that maybe you'd want a step-by-step work through of the "telescoping" method by which the sum of such a series is found. Here it is.
 

Download TelescopingSeries.mw