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.