While trying your solution, I noticed that my maple gives me a higher order of the complex part of the solution (E-6).
That's because I had Digits set to 15 and you had Digits set to 10. You can include the command Digits:= 15 right after your restart.
Furthermore I cannot plot my solution due to the following error....
Your plot command has unbalanced quotation marks, and possibly some other bad character(s). I'm surprised that it got past the parser. I recommend that you retype (not cut and paste) the whole command in a new execution group. All of the quotes need to be double quotes ("), not single quotes ('). I also recommend that you use input mode Maple Input (like the red characters that I used) rather than 2D Input, because with Maple Input both you and I can see exactly what you typed. With 2D Input it is impossible to visually distinguish a double quote from two single quotes next to each other.
Furthermore I get an error trying to find the slope of the catenary at x=0.
Use eval(Helling, x= 0).
I was wondering why maple supplies us with two solutions for the differential equation.
There is only one solution to the differential equation proper; it is only after the boundary conditions are put in that the branching occurs. You can see that the differential equation has one rather simple solution by calling dsolve without the boundary conditions:
cosh(_C1 qT + qT x)
y(x) = ------------------- + _C2
Is it because of the square root in the system?
No. When you apply the boundary conditons to the above solution to solve for the constants of integration _C1 and _C2, you get an equation of the form cosh(x) = B. Such equations generally have two solutions (for B > 0) because horizontal lines cross cosh curves twice.
Furthermore, what is the origin of the complex part of the solution?
After the boundary conditions are applied, the solution contains an expression of the form cosh(A+ln(RootOf(Q))) where Q is a quadratic with two real roots, one positive and one negative. (The nature of the roots depends on your parameter values and boundary conditions. For the case we are discussing, one is positive and one is negative.) The ln of a negative number is complex; specifically, for x < 0, ln(x) = ln(|x|) + Pi*I.
It is important to note that this is not a complex solution; we are not ignoring any part of the true solution. The Pi*I turns real after the cosh is applied to it. It is only because of the inexact floating-pointing arithmetic that we see it at all. That's why I call them spurious imaginary parts. Below, I redid the worksheet in such a way that exact arithmetic is used up to the point that the complex parts disappear.