eq = y''[x] + Epsilon*Exp[x/3]*y[x] == 0
soln = DSolve[eq, y[x], x]
дает ответ, и тогда вы даже можете сделать
Series[y[x] /. soln[[1]], {Epsilon, 0, 2}]
что некрасиво.
Чтобы добавить начальные условия, вы просто добавляете уравнения:
soln2 = DSolve[{eq, y[0]==1, y'[0]==1}, y[x], x]//Simplify
Series[y[x] /. soln2[[1]], {Epsilon, 0, 2}]//Simplify
(где я добавил //Simplify
, чтобы заставить Mathematica преобразовать его в красивую форму.