Solving system of ODEs with extra parameter

As far as I can tell, you did everything right in your second approach. Classically, Mathematica returns lists of rules as results of solving functions, but in this case, I find this tradition rather confusing and prefer to use ParametricNDSolveValue; it returns a ParametricFunction object that, when applied to a numerical parameter, returns a list of 4 InterpolatingFunction for your 4 functions.

T[θ_] = {{T11[θ], T12[θ]}, {T21[θ], T22[θ]}};
A[θ_] = {{0, E^(-I θ)/λ}, {1/36 E^(-I θ) (9 λ + 2 (-1 + λ)^2 (6 + 6 Cos[θ] + Cos[2 θ])), 0}};
sys = {T'[θ] == T[θ].A[θ]};

Tsol = ParametricNDSolveValue[{sys, T11[0] == 1, T12[0] == 0, 
    T21[0] == 0, T22[0] == 1},
   {T11, T12, T21, T22},
   {θ, 0, 2 Pi},
   {λ}
   ];

In order to obtain the numerical values for all the solutions at θ = 2 Pi for a given parameter, say λ = 0.1, you may use Through:

Through[Tsol[0.1][2. Pi]]

{-0.545795 + 1.00532 I, -1.43497 - 7.95125*10^-7 I, -0.215035 - 2.80298*10^-8 I, -0.545795 - 1.00532 I}

In order to make that into a function, you may use

f = λ \[Function] Through[Tsol[λ][2. Pi]]

For this case, you don't even need to write out the components of your matrix function:

pf = ParametricNDSolveValue[{t'[θ] == t[θ].{{0, Exp[-I θ]/λ},
                                            {Exp[-I θ] (9 λ + 2 (λ - 1)^2
                                             (6 Cos[θ] + Cos[2 θ] + 6))/36, 0}},
                             t[0] == IdentityMatrix[2]}, t, {θ, 0, 2 π}, λ,
                             Method -> "StiffnessSwitching"];

sol = pf[(3 + 4 I)/5];

ParametricPlot[ReIm[Tr[sol[θ]]], {θ, 0, 2 π}]

ParametricPlot[ReIm[Det[sol[θ]]], {θ, 0, 2 π}]

You can even make a plot where the parameter is varying:

Plot[Re[Tr[pf[Exp[I ϕ]][2 π]]], {ϕ, 0, 2 π}]

Here is how to obtain the series expansion of the matrix components. First, using a tweaked version (solving for t[2π] instead of t) of JM's formulation:

A[θ_] = {
    {0, E^(-I θ)/λ},
    {1/36 E^(-I θ) (9 λ + 2 (-1 + λ)^2 (6 + 6 Cos[θ] + Cos[2 θ])), 0}
};

pf = ParametricNDSolveValue[
    {
    t'[θ] == t[θ].A[θ], t[0] == IdentityMatrix[2]
    },
    t[2π],
    {θ, 0, 2π},
    λ
];

Then, pf will return the matrix value at $\theta = 2 \pi$. For example:

pf[1]
pf[Exp[I Pi/3]]

{{1. + 3.62818*10^-9 I, 6.82646*10^-8 - 5.73536*10^-9 I}, {1.70661*10^-8 - 1.43384*10^-9 I, 1. + 3.62818*10^-9 I}}

{{0.985595 + 1.17074 I, -0.425572 + 0.737112 I}, {-0.788363 - 1.36549 I, 0.985595 - 1.17074 I}}

Finding the series expansion is simple:

DecimalForm[Series[pf[λ], {λ, 1, 5}], {4,4}] //TeXForm

$\left( \begin{array}{cc} 1.0000 & 6.2830 \\ 1.5710 & 1.0000 \\ \end{array} \right)+\left( \begin{array}{cc} 0.0000+0.0000 i & -6.2830+0.0000 i \\ 1.5710+0.0000 i & 0.0000+0.0000 i \\ \end{array} \right) (\lambda -1)+\left( \begin{array}{cc} 0.0000-0.7869 i & 0.6691-0.0000 i \\ 0.8799+0.0000 i & -0.0000+0.7869 i \\ \end{array} \right) (\lambda -1)^2+\left( \begin{array}{cc} 0.0000+0.7869 i & -1.3380+0.0000 i \\ 0.0000+0.0000 i & 0.0000-0.7869 i \\ \end{array} \right) (\lambda -1)^3+\left( \begin{array}{cc} -0.0152-0.4524 i & 1.8410-0.0000 i \\ -0.5708-0.0000 i & -0.0152+0.4524 i \\ \end{array} \right) (\lambda -1)^4+\left( \begin{array}{cc} 0.0305+0.1179 i & -2.1780+0.0000 i \\ 0.5708-0.0000 i & 0.0305-0.1179 i \\ \end{array} \right) (\lambda -1)^5+O\left((\lambda -1)^6\right)$