Code:

pseries = Series[x/Cos[x], {x, 0, 25}];
Subscript[b, n_] := SeriesCoefficient[
pseries, n]
qseries = Series[-2/Cos[x], {x, 0, 25}];
Subscript[b, n_] := SeriesCoefficient[
pseries, n]
Subscript[d, n_] := SeriesCoefficient[
qseries, n]
Subscript[a, 0] = 1;
Subscript[a, 1] = 0;
mycoef = Table[Subscript[a, n] =
(1/(n*(n - 1)))*
(-Sum[k*Subscript[a, k]*Subscript[
b, n - k - 1], {k, 1,
n - 1}] - Sum[Subscript[a, k]*
Subscript[d, n - k - 2],
{k, 0, n - 2}]), {n, 2, 20}];
myfun[x_] = Sum[Subscript[a, n]*x^n,
{n, 0, 20}]
series = Plot[myfun[x], {x, 0, 2},
PlotStyle -> Red]
sol = NDSolve[
{Cos[x]*Derivative[2][y][x] +
x*Derivative[1][y][x] - 2*y[x] ==
0, y[0] == 1, Derivative[1][y][0] ==
0}, y, {x, 0, 2}]
numeric = Plot[Evaluate[y[x] /. sol],
{x, 0, 2}, PlotStyle -> Blue]
Show[{series, numeric}]