You need to take the limit for x to zero if you want to use the formula for the Taylor expansion. This is not an efficient way to compue the Taylor series. A more efficient method is to substitute the Taylor expansion for e^x and then perform the division. Alternatively, you can use the power series of the geometric series.
So, we have:
x^3/[exp(x) - 1] =
x^3/[x + x^2/2 + x^3/6 + x^4/24 + ...] =
x^2/[1+x/2 + x^2/6 + x^3/24 +...]
If you use the geometric series power expansion, you use that:
1/(1+y) = 1 - y + y^2 - y^3 + ...
and substitute
y = x/2 + x^2/6 + x^3/24 + ...
in here.
If you want to do a division, you can do long division. You can also write the expansion in undetermined form:
x^2/[1+x/2 + x^2/6 + x^3/24 +...] = a x^2 + b x^3 + c x^4 + ...
(note that we can immediately tell that the series expansion will start with x^2).
If you multiply both sides by 1+x/2 + x^2/6 + x^3/24 +... and equate the coefficients of the powers of x you can solve for a, b , recursively. This procedure is equivalent to long division.
Then there exist a much faster way to do division of power series that is based on the Newton-Raphson technique of findng zeroes of nonlinear equations. This works as follows. In case of division of ordinary numbers, we need to compute 1/y where y is some given number. Then this amounts to solving the equation
1/x - y = 0
for x.
Newton-Raphson then gives you successive approximants x_n that satisfy the recursion:
x_{n+1} = x_{n} - (1/x_{n} - y)/(-1/x_{n}^2) =
2x_n - y x_{n}^2
Since this algorithm does not contain any divisions, it is a bona fide division algorithm.
Now, this recursion can also be used for power series. If y is not a number but a power series and x_0 is 1/y correct to zeroth order, then the algorithm will yield the x_n as the expansion of 1/y correct to
order 2^n-1. The number of correct expansion coefficients doubles at each iteration, while in case of long divsion you only get one coefficient per step.
