Python deconvolution giving unexpected result

Below I have plotted the signal (Lifetime decay) I am trying to deconvolve from an impulse response function, i.e. the divider (IRF). So I should just get the decay a bit sharper. Here is an example of a topic I look at that gives what I need:

Understanding scipy deconvolve

Please not for my code, I am using only the peak of the divider (IRF), not the entire array sequence as shown on the image.

I am using the following code to do that:

IRF = IRF * (max(decay)/max(IRF))
# replace 0s to avoid error message 
IRF = np.where(IRF == 0, 0.1, IRF)    
decay = np.where(decay == 0, 0.1, decay)  
# take only the quotient part of the result 
deconv = scipy.signal.deconvolve(decay, IRF)[0]
# "padding" the deconvolved signal so it has the same size as the original signal 
s = int((len(decay)-(len(deconv)))/2)  ## difference on each side 
deconv_res = np.zeros(len(decay))   
end = int(len(decay)-s-1)  # final index
deconv_res[s:end] = deconv
deconv = deconv_res
# convolved normalized to decay height for plotting 
deconv_n = deconv * (max(decay)/max(deconv))   

The IRF is an array of float64, the signal is an array of uint16.

I admit I'm not so familiar with the maths of deconvolution, so I am trying blindly different things, like trying different divider functions, but nothing is producing anywhere near as expected. The last result I got looks like this (see plot of the original signal and what the signal it tried to deconvolve..)

Could anyone give me some idea if it's something in scipy.deconvolve I don't understand, what the reason could be for this strange behaviour, or even some high-level reading material that might help me out? Or if you think this problem is more physics-y than coding-related, a better place to ask such a question?

The problem is that deconvolve is a sort of polynomial division, it decomposes the output signal in $conv(h, x) + r$, if your signal is noisy it may give strange results. Also if the first sample in the inpulse response is small it tends to produce the exponentially growing output.

What I would do for this problem is the division of FFTs.

N = 2**(ceil(log2(len(IRF) + len(decay)))
filtered = ifft(fft(decay, N) / fft(IRF, N))