Electronic – Theory question about “j” imaginary unit (AC circuit analysis)

"notice that Xl is not written as jXl"

You've answered your own question. The problem is that \$X_L\$ is not written as \$j X_L\$. :-)

While it also has the dimension Ω, just like impedance, \$X\$ is the symbol for reactance, not impedance, that would be \$Z\$.

\$ Z_L = R + X_L = R + j \omega L\$

The factor \$j\$ indicates a 90° phase shift with respect to the resistance which does not have a phase shift. For an inductance the phase shift is positive, for a capacitance it's negative, since

\$ Z_C = R + X_C = R + \dfrac{1}{j \omega C} = R - j \dfrac{1}{\omega C}\$

It's because of the 90° angle between resistance and reactance that you can't add them arithmetically, like 5 Ω + 12 Ω = 17 Ω. \$X\$ and \$R\$ create a right angled triangle, so you have to apply Pythagoras to find the magnitude of the impedance:

\$ |Z_L| = \sqrt{|R|^2 + |X_L|^2} \$

so for \$R\$ = 5 Ω and \$X_L\$ = \$j\$ 12 Ω we get a magnitude of 13 Ω instead of 17 Ω.

Right, back to your calculations. The reason it seems to go in the right direction, even without \$j\$ is that you applied Pythagoras implicitely in the first line already; you wrote

\$ V_S^2 = (V_1 + V_R)^2 + V_L^2 \$

and not

\$ V_S^2 = (V_1 + V_R + V_L)^2 \$

or simply

\$ |V_S| = |V_1 + V_R + V_L| \$

That's why you won't see \$j\$ anymore; you got rid of it at the start.

The solution to a differential equation involves arbitrary constants (\$c_1\$ and \$c_2\$ in this instance) whose value cannot be determined unless the initial conditions are known. It is not necessary that one of \$c_1\$ and \$c_2\$ be zero. Both could be nonzero, and you will still observe a (single) sinusoid of the same frequency. As @pjc50 cogently observed, it is a question of phase (that is, the initial phase value). Note that if both \$a\$ and \$b\$ are nonzero, then the sum \$a\cos(\omega t) + b\sin(\omega t)\$ is not the sum of two sinusoids of different frequencies, but rather a single sinusoid of frequency \$\omega\$ rad/s of a different initial phase:

$$\begin{align*} a\cos(\omega t) + b\sin(\omega t) &= \sqrt{a^2+b^2}\left[\frac{a}{\sqrt{a^2+b^2}} \cos(\omega t) + \frac{b}{\sqrt{a^2+b^2}}\sin(\omega t)\right]\\ &= \left.\left.\sqrt{a^2+b^2}\right[\sin(\theta) \cos(\omega t) + \cos(\theta)\sin(\omega t)\right]\\ &= \sqrt{a^2+b^2}\sin(\omega t + \theta)\\ &= \sqrt{a^2+b^2}\cos(\omega t + (\theta - \pi/2)) \end{align*}$$ where \$\theta = \arctan\left(\frac{a}{b}\right)\$ and we have used the identity \$\sin(x) = \cos(x-\pi/2)\$ in getting to that last line.