Linear differential system of Bessel equations I have the following system:
\begin{cases}
g_{1}^{\prime\prime}+\rho^{-1}g_{1}^{\prime}-(1+\rho^{-2})g_{1}-f_{3} &amp; =0\\
f_{3}^{\prime\prime}+\rho^{-1}f_{3}^{\prime}-(1+9\rho^{-2})f_{3}-g_{1} &amp; =0
\end{cases}

which I rewrote as

$$\left(\begin{array}{cc}
\mathcal{L}_{1}^{(\tilde{\kappa})} &amp; -1\\
-1 &amp; \mathcal{L}_{3}^{(\tilde{\kappa})}
\end{array}\right)\left(\begin{array}{c}g_{1}\\
f_{3}\end{array}\right)=\left(\begin{array}{c}0\\0
\end{array}\right)$$

where $\mathcal L^{(\kappa)}_n$ is the Bessel operator with solutions $K_n(\kappa r)$, $I_n(\kappa r)$.

But I am not sure how to solve it.
I assume $g_3$ and $f_1$ depend on $\rho$, since you write the system in terms of Bessel operators. Combining the two equations, it can be concisely written as
\begin{equation}
\left( \mathcal{L}_3 \mathcal{L}_1 -1 \right) g_1 = 0, \tag{1}
\end{equation}
where $\mathcal{L}_\nu$ is the operator giving the modified (!) Bessel equation
\begin{equation}
 \mathcal{L}_\nu \phi = \rho^2 \phi_{\rho\rho} + \rho \phi_\rho - (\rho^2 + \nu^2)\phi = 0.
\end{equation}
I tried to find two second order Sturm-Liouville operators $D_i$ such that $(1)$ can be written in the form
\begin{equation}
 D_1 D_2 \phi = 0,
\end{equation}
such that you could infer that $D_2\phi \in \text{ker} D_1$, which would allow you to solve the inhomogeneous second order Sturm-Liouville ODE
\begin{equation}
 D_2 \phi = \psi,\quad \psi \in \text{ker }D_1
\end{equation}
by variation of parameters. Sadly, I didn't succeed in finding such Sturm-Liouville operators, but maybe you (or someone else) will have more luck.

A general idea is to use the Fourier-Bessel series to express a general solution as a sum of Bessel functions (actually, modified Bessel functions in your case, but the idea is the same). Because these functions are eigenfunctions of the differential operators $\mathcal{L}_\nu$, equation $(1)$ might give you a particularly simple relation between the series coefficients.
Solving by large computer algebra system (Mathematica 10.3.1), does not turn up Bessel functions.  The form is similar to that of liniear fourht order equations, suggesting that the solution technique is to solve, for example, your second equation for $g_1$, then substitute that into the first equation, yielding a fourth order linear equation.

[The original solution without cleanup is preserved far below.]

Solution from Mathematica is below.  The four constants of integration are $c_1$, $c_2$, $c_3$, and $c_4$.  The code to generate this output is

However, somewhat more readable is 

FullSimplify[FullSimplify[
 {f3[x], g1[x]} /. DSolve[{
    g1''[x] + ρ^-1 g1'[x] - (1 + ρ^-2) g1[x] - f3[x] == 0,
    f3''[x] + ρ^-1 f3'[x] - (1 + 9 ρ^-2) f3[x] - g1[x] == 0
    },
   {f3, g1},
   x]
 ] //. {Sqrt[ρ^4 (16 + ρ^4)] -&gt; R1, 1/Sqrt[ρ^4 (16 + ρ^4)] -&gt; 1/R1,
Sqrt[-4 R1 + 21 ρ^2 + 4 ρ^4] -&gt; R2, 1/Sqrt[-4 R1 + 21 ρ^2 + 4 ρ^4] -&gt; 1/R2, 
Sqrt[4 R1 + 21 ρ^2 + 4 ρ^4] -&gt; R3,  1/Sqrt[4 R1 + 21 ρ^2 + 4 ρ^4] -&gt; 1/R3}
] // TraditionalForm // TeXForm


giving (where I have made one manual change, moving the denominator to a fraction in front) $\{f_3(x), g_1(x)\} = {}$

$\left\{\frac{1}{4 R_1 R_2 R_3} \left( R_3 \left(R_1 \left(\rho  \left(2 c_2 \rho +c_1\right)+c_1
   R_2\right)-\rho ^2 \left(\rho  \left(\rho  \left(\rho  \left(2 c_4 \rho
   +c_3\right)+8 c_2\right)+4 c_1\right)+R_2 \left(c_3 \rho ^2+4
   c_1\right)\right)\right) e^{\frac{x \left(R_2-\rho \right)}{2 \rho
   ^2}}+R_3 \left(\rho ^3 \left(\rho  \left(\rho  \left(2 c_4 \rho
   +c_3\right)+8 c_2\right)+4 c_1\right)+R_2 \left(-c_3 \rho ^4-4 c_1 \rho
   ^2+c_1 R_1\right)-\rho  R_1 \left(2 c_2 \rho +c_1\right)\right)
   e^{-\frac{x \left(\rho +R_2\right)}{2 \rho ^2}}+R_2 \left(\rho ^3
   \left(-\left(\rho  \left(\rho  \left(2 c_4 \rho +c_3\right)+8
   c_2\right)+4 c_1\right)\right)+R_3 \left(c_3 \rho ^4+4 c_1 \rho ^2+c_1
   R_1\right)-\rho  R_1 \left(2 c_2 \rho +c_1\right)\right) e^{-\frac{x
   \left(\rho +R_3\right)}{2 \rho ^2}}+R_2 \left(\rho ^3 \left(\rho 
   \left(\rho  \left(2 c_4 \rho +c_3\right)+8 c_2\right)+4 c_1\right)+R_3
   \left(c_3 \rho ^4+4 c_1 \rho ^2+c_1 R_1\right)+\rho  R_1 \left(2 c_2 \rho
   +c_1\right)\right) e^{\frac{x \left(R_3-\rho \right)}{2 \rho ^2}}\right),\frac{1}{4 R_1 R_2 R_3} \left( R_3 \left(\rho ^3 \left(\rho  \left(\rho  \left(2 c_2 \rho
   +c_1\right)-8 c_4\right)-4 c_3\right)+R_2 \left(-c_1 \rho ^4+4 c_3 \rho
   ^2+c_3 R_1\right)-\rho  R_1 \left(2 c_4 \rho +c_3\right)\right)
   e^{-\frac{x \left(\rho +R_2\right)}{2 \rho ^2}}+R_3 \left(\rho ^3
   \left(\rho  \left(8 c_4-\rho  \left(2 c_2 \rho +c_1\right)\right)+4
   c_3\right)+R_2 \left(-c_1 \rho ^4+4 c_3 \rho ^2+c_3 R_1\right)+\rho  R_1
   \left(2 c_4 \rho +c_3\right)\right) e^{\frac{x \left(R_2-\rho \right)}{2
   \rho ^2}}+R_2 \left(\rho ^3 \left(\rho  \left(8 c_4-\rho  \left(2 c_2
   \rho +c_1\right)\right)+4 c_3\right)+R_3 \left(c_1 \rho ^4-4 c_3 \rho
   ^2+c_3 R_1\right)-\rho  R_1 \left(2 c_4 \rho +c_3\right)\right)
   e^{-\frac{x \left(\rho +R_3\right)}{2 \rho ^2}}+R_2 \left(\rho ^3
   \left(\rho  \left(\rho  \left(2 c_2 \rho +c_1\right)-8 c_4\right)-4
   c_3\right)+R_3 \left(c_1 \rho ^4-4 c_3 \rho ^2+c_3 R_1\right)+\rho  R_1
   \left(2 c_4 \rho +c_3\right)\right) e^{\frac{x \left(R_3-\rho \right)}{2
   \rho ^2}} \right)\right\}$

where

$$\begin{align}
    R_1 &amp;= \sqrt{\rho^4(16+\rho^4)} \\
    R_2 &amp;= \sqrt{4 \rho^4 + 21 \rho^2 - 4 R_1} \\
    R_3 &amp;= \sqrt{4 \rho^4 + 21 \rho^2 + 4 R_1}
\end{align}$$

Note that $x$ appears only in the exponents.  The method of undetermined coefficients applied to these exponentials (and "$1$") should clean this up substantially.  However, my original observation is unchanged: these exponents are from the characteristic equation arising from the fourth order equations you get by solving your first equation for $f_3$ and substituting into your second equation and vice versa for $g_1$.

Original mess below here.  Note that there is some room for compaction by extracting frequently repeated (constant) subexpressions.  The four constants of integration are $c_1$, $c_2$, $c_3$, and $c_4$.  The code to generate this output is

FullSimplify[
 {f3[x], g1[x]} /. DSolve[{
    g1''[x] + ρ^-1 g1'[x] - (1 + ρ^-2) g1[x] - f3[x] == 0,
    f3''[x] + ρ^-1 f3'[x] - (1 + 9 ρ^-2) f3[x] - g1[x] == 0
    },
   {f3, g1},
   x]
 ]
% // TeXForm


$\{ f_3(x), g_1(x)\} = \\
\left\{\frac{e^{\frac{x \left(\sqrt{4 \rho ^4+21 \rho ^2-4
   \sqrt{\rho ^4 \left(\rho ^4+16\right)}}-\rho \right)}{2 \rho ^2}} \sqrt{4
   \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(-2 c_4
   \rho ^6-c_3 \rho ^5-\left(8 c_2+\sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho
   ^4 \left(\rho ^4+16\right)}} c_3\right) \rho ^4-4 c_1 \rho ^3+\left(2
   \sqrt{\rho ^4 \left(\rho ^4+16\right)} c_2-4 \sqrt{4 \rho ^4+21 \rho ^2-4
   \sqrt{\rho ^4 \left(\rho ^4+16\right)}} c_1\right) \rho ^2+\sqrt{\rho ^4
   \left(\rho ^4+16\right)} c_1 \rho +\sqrt{\rho ^4 \left(\rho ^4+16\right)}
   \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   c_1\right)+e^{-\frac{x \left(\rho +\sqrt{4 \rho ^4+21 \rho ^2-4
   \sqrt{\rho ^4 \left(\rho ^4+16\right)}}\right)}{2 \rho ^2}} \sqrt{4 \rho
   ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(2 c_4 \rho
   ^6+c_3 \rho ^5+\left(8 c_2-\sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}} c_3\right) \rho ^4+4 c_1 \rho ^3-2 \left(2
   \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   c_1+\sqrt{\rho ^4 \left(\rho ^4+16\right)} c_2\right) \rho ^2-\sqrt{\rho
   ^4 \left(\rho ^4+16\right)} c_1 \rho +\sqrt{\rho ^4 \left(\rho
   ^4+16\right)} \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}} c_1\right)-e^{-\frac{x \left(\rho +\sqrt{4 \rho ^4+21 \rho
   ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}\right)}{2 \rho ^2}} \sqrt{4
   \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(2 c_4
   \rho ^6+c_3 \rho ^5+\left(8 c_2-\sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho
   ^4 \left(\rho ^4+16\right)}} c_3\right) \rho ^4+4 c_1 \rho ^3+\left(2
   \sqrt{\rho ^4 \left(\rho ^4+16\right)} c_2-4 \sqrt{4 \rho ^4+21 \rho ^2+4
   \sqrt{\rho ^4 \left(\rho ^4+16\right)}} c_1\right) \rho ^2+\sqrt{\rho ^4
   \left(\rho ^4+16\right)} c_1 \rho -\sqrt{\rho ^4 \left(\rho ^4+16\right)}
   \sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   c_1\right)+e^{\frac{x \left(\sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}}-\rho \right)}{2 \rho ^2}} \sqrt{4 \rho ^4+21
   \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(2 c_4 \rho ^6+c_3
   \rho ^5+\left(8 c_2+\sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}} c_3\right) \rho ^4+4 c_1 \rho ^3+2 \left(2 \sqrt{4 \rho
   ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} c_1+\sqrt{\rho ^4
   \left(\rho ^4+16\right)} c_2\right) \rho ^2+\sqrt{\rho ^4 \left(\rho
   ^4+16\right)} c_1 \rho +\sqrt{\rho ^4 \left(\rho ^4+16\right)} \sqrt{4
   \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   c_1\right)}{4 \sqrt{\rho ^4 \left(\rho ^4+16\right)} \sqrt{4 \rho ^4+21
   \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \sqrt{4 \rho ^4+21 \rho
   ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}},\frac{e^{-\frac{x
   \left(\rho +\sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}}\right)}{2 \rho ^2}} \sqrt{4 \rho ^4+21 \rho ^2+4
   \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(2 c_2 \rho ^6+c_1 \rho
   ^5-\left(\sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}} c_1+8 c_4\right) \rho ^4-4 c_3 \rho ^3+\left(4 \sqrt{4
   \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} c_3-2
   \sqrt{\rho ^4 \left(\rho ^4+16\right)} c_4\right) \rho ^2-\sqrt{\rho ^4
   \left(\rho ^4+16\right)} c_3 \rho +\sqrt{\rho ^4 \left(\rho ^4+16\right)}
   \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   c_3\right)+e^{\frac{x \left(\sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}}-\rho \right)}{2 \rho ^2}} \sqrt{4 \rho ^4+21
   \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(-2 c_2 \rho
   ^6-c_1 \rho ^5+\left(8 c_4-\sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}} c_1\right) \rho ^4+4 c_3 \rho ^3+2 \left(2
   \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   c_3+\sqrt{\rho ^4 \left(\rho ^4+16\right)} c_4\right) \rho ^2+\sqrt{\rho
   ^4 \left(\rho ^4+16\right)} c_3 \rho +\sqrt{\rho ^4 \left(\rho
   ^4+16\right)} \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}} c_3\right)+e^{-\frac{x \left(\rho +\sqrt{4 \rho ^4+21 \rho
   ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}\right)}{2 \rho ^2}} \sqrt{4
   \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} \left(-2 c_2
   \rho ^6-c_1 \rho ^5+\left(\sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}} c_1+8 c_4\right) \rho ^4+4 c_3 \rho ^3-2
   \left(2 \sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}} c_3+\sqrt{\rho ^4 \left(\rho ^4+16\right)} c_4\right) \rho
   ^2-\sqrt{\rho ^4 \left(\rho ^4+16\right)} c_3 \rho +\sqrt{\rho ^4
   \left(\rho ^4+16\right)} \sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}} c_3\right)+e^{\frac{x \left(\sqrt{4 \rho ^4+21
   \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}-\rho \right)}{2 \rho
   ^2}} \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}}
   \left(2 c_2 \rho ^6+c_1 \rho ^5+\left(\sqrt{4 \rho ^4+21 \rho ^2+4
   \sqrt{\rho ^4 \left(\rho ^4+16\right)}} c_1-8 c_4\right) \rho ^4-4 c_3
   \rho ^3+\left(2 \sqrt{\rho ^4 \left(\rho ^4+16\right)} c_4-4 \sqrt{4 \rho
   ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho ^4+16\right)}} c_3\right) \rho
   ^2+\sqrt{\rho ^4 \left(\rho ^4+16\right)} c_3 \rho +\sqrt{\rho ^4
   \left(\rho ^4+16\right)} \sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4
   \left(\rho ^4+16\right)}} c_3\right)}{4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)} \sqrt{4 \rho ^4+21 \rho ^2-4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}} \sqrt{4 \rho ^4+21 \rho ^2+4 \sqrt{\rho ^4 \left(\rho
   ^4+16\right)}}}\right\}$