I don't really get a solution to a problem:
Problem: $y''+4y=\frac{1}{\sin(2x)}$
We calculate the general homogeneouse solution using the ansatz $y=e^{\lambda x}$
So we get:
$$z^2+4=(z+2i)(z-2i)$$
The solution is $y(x)=A\cos(2x)+B\sin(2x)$ with $A,B \in \mathbb R$
Let $A,B$ be two functions and
$$y_0(x)=A(x)\cos(2x)+B(x)\sin(2x)$$
such that
$$\begin{cases}A' (x)\cos(2x)+B'(x)\sin(2x)=0\\-2A'(x)\sin(2x)+2B'\cos(2x)=\frac{1}{\sin(2x)}\end{cases} \tag{1}$$
For all $x\in \mathbb R$, the matrix
$$R(x)=\begin{pmatrix}\cos(2x)&\sin(2x)\\-2\sin(2x)&2\cos(2x)\end{pmatrix}\tag{2}$$
Then because of (1) and (2) we get:
$$\begin{cases}A'(x)=\frac{-1}{2}\\B'(x)=\frac{\cos(2x)}{2\sin(2x)}\end{cases}$$
is invertible. So:
$$R(x)^{-1}=\frac{1}{2}\begin{pmatrix}2\cos(2x)&-2\sin(2x)\\ 2\sin(2x)&\cos(2x)\end{pmatrix}$$
Because of (1) and (2) we get:
$$\begin{cases}A'(x)=\frac{-1}{2}\\ B'(x)=\frac{\cos(2x)}{2\sin(2x)}\end{cases}$$
and
$$\begin{cases}A(x)=\frac{-x}{2}\\ B(x)=\frac{1}{4}\log|\sin(2x)|\end{cases}$$
so we get
$$y(x)=A\cos(2x)+B\sin(2x)-\frac{x}{2}\cos(2x)+\frac{1}{4}\log|\sin(2x)|\sin(2x)$$
Question: I do think I get the general idea of variation of constants. Basically, we are variating the constants to make our homogeneous solution space fit the inhomogenity. But I absolutely don't see what was done at (1). How did they get that?
$$y''+4y=\frac{1}{\sin(2x)}$$ You found the general solution of the homogeneous ODE : $$y_h''+4y_h=0 \qquad;\qquad y_h=A\cos(2x)+B\sin(2x)$$ $A$ and $B$ are constants.
Then you replace the constants by unknown functions $A(x)$ and $B(x)$ and you look for a particular solution $y_p$ of the non-homogeneous ODE , thus satisfying : $$\begin{cases} y_p''+ 4y_p=\frac{1}{\sin(2x)} \\ y_p=A(x)\cos(2x)+B(x)\sin(2x) \end{cases}$$ $\left(A(x)\cos(2x)+B(x)\sin(2x)\right)''+4\left(A(x)\cos(2x)+B(x)\sin(2x)\right)=\frac{1}{\sin(2x)} $
Note that you have two unknows $A(x)$ and $B(x)$ but only one equation. So, you are allowed to add an arbitrary relationship between $A(x)$ and $B(x)$. This will be useful later for simplification.
$y_p'=A'\cos(2x)+B'\sin(2x)-2A\sin(2x)+2B\cos(2x)$
$y_p''=A''\cos(2x)+B''\sin(2x)-4A'\sin(2x)+4B'\cos(2x)-4A\cos(2x)+4B\sin(2x)$
$$y_p''+ 4y_p=A''\cos(2x)+B''\sin(2x)-4A'\sin(2x)+4B'\cos(2x)=\frac{1}{\sin(2x)}$$
Using
$\left(A'\cos(2x)+B'\sin(2x)\right)'=A''\cos(2x)+B''\sin(2x)-2A'\sin(2x)+2B'\cos(2x)$
the equation becomes :
$$\left(A'\cos(2x)+B'\sin(2x)\right)'-2A'\sin(2x)+2B'\cos(2x)=\frac{1}{\sin(2x)}$$
At this point, the advantage to introduce an relationship between $A(x)$ and $B(x)$ appears :
Let : $$A'\cos(2x)+B'\sin(2x)=0$$
The equation becomes much simpler : $$-2A'\sin(2x)+2B'\cos(2x)=\frac{1}{\sin(2x)}$$
So, you get the system $(1)$ of two ODEs allowing to compute $A(x)$ and $B(x)$ and then the particular solution $y_p(x)$ and the general solution $ y(x)=y_h(x)+y_p(x)$.