If you were to use an if-statement, what would be the condition---the expression that's true or false---that would start the statement?

You're overthinking it, there's no such condition, and your procedure only needs one statement, based on this hint: If f is a function, and a a value in its domain, and the coordinate axes are named x and y, then the tangent line is y = D(f)(a)*(x-a)+f(a).

Anywhere in Maple, if a name is enclosed by single forward quotes (aka aposthropes), then it'll be treated as a name even if it has been assigned a value. So,

TangentLine('l', point('A', 0, 1), circle('c', [point('C', 0, 1), 1]));

When quotes are used this way, they're called unevaluation quotes.

In Maple 2019, the semantics have changed thus: Certain bound variables that were implicitly not local in prior versions are now implicitly local (unless explicitly made otherwise). And the syntax change is that (as shown by VV) it's been made easier to declare things local.

While I don't believe in heeding every warning issued by Maple, I think that it's a cruddy programming practice to ignore this one. Declaring all variables will help immensely with your debugging if you ever misspell a variable's name. I realized that all the way back in my Fortran days (1979). That compiler let you leave all integer and real variables implicit (based on the first letter of the name) if you so chose.

I agree with all that Acer wrote in this thread. But there's an additional tidbit of information that he didn't mention that I think greatly simplifies determining whether an argument was passed to a parameter:

• The default value assigned to an optional parameter (whether positional or keyword) need not be of the type declared to be  allowed for arguments to that parameter.

So, if you use such a "special" default value, and the parameter has that value, then it's guaranteed that no argument was passed to it.

For a.e. x and t, there's nothing to suggest that the limit of the terms is 0. Indeed, the terms appear to be a dense subset of -1..1. If these appearances hold true, then the series fails the most basic test of convergence, the "divergence test".

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.