If you leave a space when you type sec (theta) then Maple thinks you are trying to multiply (because juxtaposition with a space is an implicit multiplication---yes I know lots of people wanted that "feature" but it causes lots of problems).

Anyway.  If you have r = some expression in theta, which you can get by saying (notice, no r(theta), although you *can* do that; I don't recommend it, though).

r := sec( theta-Pi/2*(2/Pi *theta+Pi/4) ) 

which (as was pointed out) is just a constant function r = sec( Pi^2/8 ) which is sort of weird but ok if that's what you want, 

then you can plot the cartesian version of this by using the definition: x = r*cos(theta) and y = r*sin(theta).  This is most easily done in Maple by using the parametric plots feature (which, IMO, is under-utilized in teaching; it's very useful, very expressive, and reinforces a bunch of important concepts).  You could even animate the parameterization.

plot( [r*cos(theta), r*sin(theta), theta=-Pi..Pi] )

which, in this case, produces a circle with radius r=sec(Pi^2/8).  This may or may not have been what you meant.

If you have an actual r as a function of theta, maybe 

r := sec( theta - Pi/4 )

then the same thing works. Screenshot (showing Maple gratuitously changes the sec to a csc, but fine).

The use of "solve" will first try an analytic solution (...sometimes) and this can be expensive and may have numerical difficulties---the analytic formulas for solving a cubic and quartic can be numerically unstable (so can the quadratic formula).

So if you know you want a numerical solution for a polynomial, then use fsolve.

Separation into real and imaginary parts is straightforward, as is plotting the roots.

fsolvedemo.mw

Dan's work is good, as is John Carlos Baez's on the topic----but Maple can play here too.  Working with matrices and polynomials together is possible.  See my Maple Transactions article on the topic (mapletransactions.org) https://doi.org/10.5206/mt.v1i1.14039

(that's actually a workbook, which has quite a bit of code in it that you can alter at will)

For more fun you can see the links at bohemianmatrices.com as well.

There's a lot to learn here.  Welcome to the Bohemian world!

If you graph U(t)/U_0 versus time, you will see the same shape of graph no matter what U_0 turns out to be.  We say that the function has been "nondimensionalized" because while U(t) and U_0 are measured in volts (or whatever) their ratio is a pure number.

Such graphs are very common in engineering and physics.  Indeed, if you want to nondimentionalize the frequency as well, you can do that: put tau = omega*t (in your example omega is 2) and now you have only one parameter left, namely the damping coefficient: you have a numerical value for that, 1/10, but you can replace t/10 with delta*tau = tau/20 because omega is 2, and now plot U/U0 versus tau AND delta in a 3d plot and see the effect of various damping coefficients all at once in one graph.

With explicit damping and frequency as a one-dimensional graph (t in, say, seconds) (Screenshots here, link to the worksheet below)

and a 3d graph where all axes represent pure numbers: U/U0 versus tau and delta:

Nondimensionalize.mw

HarmonicOscillator.mw

The above worksheet contains an explicit solution to the harmonic oscillator equation y'' + sin(y) = 0.  In theory, Maple should be able to turn its integrals into the explicit solution shown there (which use Jacobian elliptic functions).  

But it is nice to know that Maple can at least verify numerically that this analytical solution is a solution.

We all have our preferred methods for this.  Mine is to use eigenvalue techniques (I admit that for large polynomials there are better methods---I should mention MPSolve).  If your polynomial actually is the characteristic polynomial of a matrix and you know the matrix already, then use that, and LinearAlgebra:-Eigenvalues on the floating-point version of your matrix.

If you do not know a matrix whose eigenvalues are the roots of the polynomial, then you can select from an infinite number of "companion" matrices.  See the attached worksheet.

PolynomialZeroFindingByEigenvalues.mw

The equations you give are an example of a nonlinear system of univariate ordinary differential equations; more, it is a boundary value problem, not just an initial value problem.  I didn't check if you included constraints in your equations, making it a Differential Algebraic equation, which is even harder.  It is very rare that nonlinear BVP for ODE can be solved symbolically ("dsolve" is working away at it as I type this, just in case a miracle occurs).  There are good numerical techniques for solving such, however; given boundary conditions and numerical values for the parameters.  See ?dsolve, bvp for instructions.

Since you haven't included many details of your problem, I'm going to take the opportunity to give a generic answer.  Simplification is actually very hard; in fact, for most classes of expressions, provably impossible: one cannot write a computer program that will reliably recognize when the expressions simplify to zero.  This is a classic result due to Dan Richardson.

Fine, we still want to simplify specific expressions.  What can we do?  Sometimes factor is good, but sometimes not (consider x^n - 1 for very large n, which has many factors).  Sometimes expand is good, but sometimes not (consider (x-1)^n, for very large n).  Sometimes collecting terms is good; one of my favourite techniques is to collect an expression in one variable, and apply a function (frequently "factor") to each coefficient: collect( bigmess, z, factor ) will give you a polynomial in z with factored coefficients (if bigmess is a polynomial in z).

One useful technique for series solutions of PDE is the use of Veil (from the LargeExpressions package). From the help page for Veil:

Which technique to use on a given problem is, of course, the hard question.  Doing the right things in the right order is an art, and really you have to develop the expertise yourself for the particular class of problem you are interested in.

Try amputating your problem, maybe, to make a smaller simpler problem to experiment on (e.g. set a bunch of parameters to zero, just so you can play around).

If you expand first, to get cos(3)*cos(sin(x)) - sin(3)*sin(sin(x)), then because there's a definite integral for Bessel functions you can find the integral from 0 to Pi (not the integral from 0 to 5) as cos(3)*Pi*BesselJ(0,1) - sin(3)*Pi*StruveH(0,1).  I do not know if there is an expression for the definite integral to another limit than Pi (and neither, apparently, does Maple).  This may not help you, but I find it interesting.