The corresponding homogeneous equation
[tex]\dfrac{\mathrm d^2y}{\mathrm dx^2}+y=0[/tex]
has characteristic equation
[tex]r^2+1=0[/tex]
which admits two roots, [tex]r=\pm i[/tex], so the characteristic solution is
[tex]y_c=C_1\cos x+C_2\sin x[/tex]
So we know two fundamental solutions, [tex]y_1=\cos x[/tex] and [tex]y_2=\sin x[/tex]. These two solutions have Wronskian determinant
[tex]W(y_1,y_2)=\begin{vmatrix}\cos x&\sin x\\-\sin x&\cos x\end{vmatrix}=\cos^2x+\sin^2x=1[/tex]
so they are linearly independent.
We use variation of parameters to find a particular solution [tex]y_p=u_1y_1+u_2y_2[/tex], where
[tex]u_1=-\displaystyle\int\sin x(\tan x+3x-1)\,\mathrm dx=-\ln|\sec x+\tan x|-2\sin x+(3x-1)\cos x[/tex]
[tex]u_2=\displaystyle\int\cos x(\tan x+3x-1)\,\mathrm dx=2\cos x+(3x-1)\sin x[/tex]
[tex]\implies y_p=-\cos x\ln|\sec x+\tan x|-2\sin x\cos x+(3x-1)\cos^2x+2\sin x\cos x+(3x-1)\sin^2x[/tex]
[tex]\implies y_p=-\cos x\ln|\sec x+\tan x|+3x-1[/tex]
[tex]\implies\boxed{y=C_1\cos x+C_2\sin x-\cos x\ln|\sec x+\tan x|+3x-1}[/tex]
