Solution (source code)

= Solution

Multiply the equation for $y_i$ by $y_j$ and the equation for $y_j$ by $y_i$, then subtract. The $q$ terms cancel and the <derivative> terms combine into a <Wronskian> boundary <derivative>:
$$
\frac d{dx}[p(y_jy_i'-y_iy_j')]+(\lambda_i-\lambda_j)wy_iy_j=0.
$$
Integrating over the interval eliminates the boundary term because both <eigenfunctions> vanish at both endpoints. Since their <eigenvalues> differ,
$$
\boxed{\int_0^1w(x)y_i(x)y_j(x)\,dx=0.}
$$
This is <orthogonality> in a <weighted inner product> in the regular Sturm-Liouville setting; the endpoint values must justify the <integration by parts>.