Is this a valid partial fraction decomposition?

Write $\dfrac{4x+1}{x^2 - x - 2}$ using partial fractions.

$$ \frac{4x+1}{x^2 - x - 2} = \frac{4x+1}{(x+1)(x+2)} = \frac{A}{x+1} + \frac{B}{x-2} = \frac{A(x-2)+B(x+1)}{(x+1)(x-2)}$$

$$4x+1 = A(x-2)+B(x+1)$$

$$x=2 \Rightarrow 4 \cdot2 + 1 = A(0) + B(3) \Rightarrow B = 3$$

$$x = -1 \Rightarrow 4(-1) +1 = A(-3)+ B(0) \Rightarrow A = 1$$

Thus,

$$\frac{4x+1}{x^2-x-2} = \frac{1}{x+1} + \frac{3}{x-2}\textrm{.}$$

The substitution of $x$ ($x = 2, -1$) is a common method to find out the coefficient of the partial fractions. However, the equation on the third line is obtained by multiplying $(x+1)(x-2)$, which is assumed to be nonzero. Here we have a contradiction. Furthermore, the original function is not defined at $x=-1,2$.

How can we substitute these value for $x$? So is this method valid and rigorous? How to modify it so that it is rigorous?

This is a classical source of confusion. But the the problem is avoided completely if one takes the limit of both sides as $x \to 2$, instead of substituting $x=2$.

Clearly the result is the same, but with this formulation there is no evaluation of functions at points where they are undefined.

$A(x-2)+B(x+1) = 4x+1$ gives us a system of two equations:

$$\begin{split} A+B &= 4 \\ -2A + B &= 1 \end{split}$$

But really the first equation as really $(A+B)x = 4x$, so the substitution $x=2$ simply gives us the system:

$$\begin{split} 2A+2B &= 8 \\ -2A + B &= 1 \\ \end{split}$$

Which conveniently sums to $3B = 9$, $B=3$. It's just a choice of substitution that makes it easy to solve the system. Any such substitution is valid, just that $x=2$ and $x=-1$ happen to be particularly good choices.