<s>[INST] Is there any other approach available without trig substitution? 
Let $f\colon [0, 1] \to \mathbb{R}$ be a continuous function such that for any $x, y \in [0, 1]$, $xf(y) + yf(x) \le  1$. Find
maximum value of $\int_0^1f(x)\,dx$.

Method given was using $x= \sin\theta$, we get $$\int_0^1 f(x)\, dx = \int_0^{\pi/2} f(\sin\theta)\cos\theta\, d\theta \tag{I}$$and similarly by $x= \cos\theta$, we would get $$\int_0^1 f(x)\,dx = \int_0^{\pi/2} f(\cos\theta)\sin\theta \,d\theta \tag{II},$$ adding $(\text{I})$ and $(\text{II})$ we get max value of the required integral to be $\pi/4$. Any other approach because the idea is very tricky and one may not realize for this type of substitution.
 [/INST] I don't have an alternative approach to suggest. The first thing that occurred to me, which was to integrate over the two-dimensional region $[0,1]^2$, doesn't work:
\begin{align}
1=\int_0^1\int_0^1 1\,dx\,dy\ge\int_0^1\int_0^1(yf(x)+xf(y))\,dx\,dy&=2\int_0^1y\,dy\int_0^1 f(x)\,dx\\
&=\left[y^2\right]_0^1\int_0^1f(x)\,dx\\
&=\int_0^1f(x)\,dx.
\end{align}
My upper bound is not as tight as yours: $1>\pi/4\approx0.785$.
What you have done is to cleverly bound the integral by restricting the integration to a single contour, namely the arc of the unit circle that lies in the first quadrant. Your bound is, in fact, the best possible, as the example $f(x)=\sqrt{1-x^2}$ shows. First confirm that your maximum is attained:
\begin{align}
\int_0^1\sqrt{1-x^2}\,dx&=\int_0^{\pi/2}\sqrt{1-\sin^2\theta}\cos\theta\,d\theta\\
&=\int_0^{\pi/2}\left(\frac{1}{2}\cos2\theta+\frac{1}{2}\right)\,d\theta\\
&=\frac{\pi}{4}.
\end{align}
To verify that $f$ satisfies the required condition, observe that the function
$$
g(x,y)=y\sqrt{1-x^2}+x\sqrt{1-y^2}
$$
achieves a maximum of $1$ for $(x,y)\in[0,1]^2$ all along your contour. You can check this by observing that the equations
\begin{align}
0&=\frac{\partial}{\partial x}\left[
y\sqrt{1-x^2}+x\sqrt{1-y^2}\right]=\sqrt{1-y^2}-\frac{xy}{\sqrt{1-x^2}}\\
0&=\frac{\partial}{\partial y}\left[
y\sqrt{1-x^2}+x\sqrt{1-y^2}\right]=\sqrt{1-x^2}-\frac{xy}{\sqrt{1-y^2}}
\end{align}
are satisfied when $x^2+y^2=1$. To see that it is a maximum that is achieved at this locus of critical points, switch to polar coordinates,
$$
h(r,\theta)=g(r\cos\theta,r\sin\theta),
$$
and compute
\begin{align}
\frac{\partial}{\partial r}h(r,\theta)=&\left(\sin\theta\sqrt{1-r^2\cos^2\theta}+\cos\theta\sqrt{1-r^2\sin^2\theta}\right)\\
&\times\left(1-\frac{r^2\sin\theta\cos\theta}{\sqrt{1-r^2\cos^2\theta}\sqrt{1-r^2\sin^2\theta}}\right),\\
\frac{\partial}{\partial\theta}h(r,\theta)=&r\cos\theta\sqrt{1-r^2\cos^2\theta}+\frac{r^3\sin^2\theta\cos\theta}{\sqrt{1-r^2\cos^2\theta}}\\
&-r\sin\theta\sqrt{1-r^2\sin^2\theta}-\frac{r^3\cos^2\theta\sin\theta}{\sqrt{1-r^2\sin^2\theta}}.
\end{align}
One can check that, for fixed $\theta$,  $\frac{\partial h(r,\theta)}{\partial r}$ is positive when $r<1$, zero when $r=1$, and negative when $r>1$, while, for $r$ fixed to $1$, $\frac{\partial h(r,\theta)}{\partial\theta}$ is zero, as  expected.

</s>