The integrals you have to compute are not that difficult. First of all, the spatial parts of your wavefunctions are real (this is wise and always possible in one dimension), and therefore complex conjugates just do not matter. So, the remaining integral you have on the rhs is (up to some expression dependent on t (S(t)) which we will leave for now):
$$ I(t) = S(t) \int_0^a\Phi_n(x)\Phi_m(x) dx = S(t) \frac{2}{a} \int_0^a \operatorname{sin}\Big( \frac{n\pi x}{a} \Big) \operatorname{sin}\Big( \frac{m\pi x}{a} \Big)dx$$
Now, you have to make a substitution $ t = \frac{\pi x}{a} $ to get
$$ I(t) = S(t) \frac{2}{\pi} \int_0^\pi \operatorname{sin} ( n x ) \operatorname{sin} ( m x )dx$$
First of all, if either n or m are 0 then I(t) = 0 and your initial state is normalized.
Now you integrate by parts (or by mathematica if you can't by parts) to obtain
$$ I(t) = S(t) \frac{2}{\pi} [ -(\frac{1}{m} \sin{n x} \cos{m x}) |_0^\pi + \frac{n}{m} \int_0^\pi \cos{n x} \cos({m x}) dx]$$
The first term is 0 (cause sin is 0 in all multiplicities of $\pi$). Now we can proceed with integration by parts of the second term.
$$ I(t) = S(t) \frac{2}{\pi} ( \frac{n}{m} \int_0^\pi \cos{n x} \cos({m x}) dx) = S(t) \frac{2}{\pi} [ (\frac{n}{m^2} \cos{n x} \sin{m x}) |_0^\pi $$
$$ + \frac{n^2}{m^2} \int_0^\pi \sin{n x} \sin({m x}) dx] = \frac{n^2}{m^2} I(t)$$
So, there are two possibilities - either n=m or I(t)=0 (n=-m is excluded as negative labels are not in the basis because sin(-nx) = -sin(nx) - linear dependency). What you get is ORTHOGONALITY RELATION - if eigenvectors have different labels (e.g. 1 and 2 as in your case) they are orthogonal in $L^2([0,a])$ which happens to be your Hilbert space. There are theorems which state that this is always true for Hamiltonian eigenstates - corresponding theory of differental operators is known as Strum-Liouville theory (for finite dimensional Hilbert spaces this is a trivial property of self-adjoined operators).
Now, to the second part of your question. The modulus squared of wavefunction is by definition probability density of finding a particle in an interval $dx$
So, the probability of finding particle in [0,a/2] is simply
$$\int_0^{a/2} |f(x,t)|^2 dx = \frac{1}{2}\int_0^{a/2} |\Phi_1(x)|^2 + |\Phi_2(x)|^2 dx + $$
$$ \frac{1}{2}\int_0^{a/2} \Phi_1(x) \Phi_2(x) \operatorname{exp} \Big( (-1^2+2^2)\frac{-i\pi^2 h t}{2ma^2} \Big)+ \Phi_2(x) \Phi_1(x) \operatorname{exp} \Big( (-2^2+1^2) \frac{-i\pi^2 h t}{2ma^2} \Big) $$
$$ = \frac{1}{2}\int_0^{a/2} |\Phi_1(x)|^2 + |\Phi_2(x)|^2 dx + \cos{\frac{3\pi^2 h t}{2ma}} \int_0^{a/2} \Phi_1(x) \Phi_2(x) dx $$
First integral is simply 1 because we take half of the area of $sin(x)^2$ and $sin(2x)$, and add it. Second integral can be easily computed using double 'by parts' formula from above and replacing $\pi$ with $\pi/2$ in integral limits. One gets $\frac{4}{3 \pi}$
Therefore finally
$$\int_0^{a/2} |f(x,t)|^2 dx = \frac{1}{2}+\frac{4}{3 \pi} \cos{\frac{3\pi^2 h t}{2ma}} $$
As $\frac{4}{3 \pi} < \frac{1}{2}$ the result makes sense. As the body is in mixture of two states, the probability is no longer constant in time.
