Solving Volterra integral equation of first kind with a Gaussian diffusive evolution kernel I am trying to solve following Voltera integral equation for $P(t|t')$ numerically:
$$
\rho(1,t|0,t') = \int_{t'}^{t} dt'' \rho(1,t|1,t'') P(t''|t')
$$
where 
$$
\rho(x,t|x',t') = \frac{1}{\sqrt{2\pi}\Sigma(t-t')}\exp(-\frac{(x-\bar{x}(t|x',t'))^2}{2\Sigma^2(t-t')}),
$$
given $\Sigma^2(t)=\frac{\sigma^2}{2\tau} [1- e^{-2t/\tau}]$ and $\bar{x}(t|x',t')$ is the solution of 
$$
(1+\frac{d}{dt})\bar{x}(t|x',t') = \mu(t).
$$
with initial condition of $\bar{x}(t|x',t)=x'$.
The problem to solve this equation numerically is that $\rho(x,t|x',t')$ is diverging while $t\rightarrow t'$ and $\rho$ becomes a delta function.
I have no idea where to begin. A solution would be appreciated, but absent that, advice on how to attack this problem would also be nice. 
Background: I believe $P$ is the first passage time of Ornstein–Uhlenbeck process form the boundary condition of $x=1$. 
Edit1: To see the problem, we may assume $t'=0$ and differentiate the main equation with respect to $t$ to obtain
$$
P(t|0) \rho(1,t|1,t) + \int_0^t \frac{\partial }{\partial t} dt'' \rho(1,t|1,t'') P(t|t'') = \frac{\partial }{\partial t} \rho(x,t|x',0)
$$
and clearly $\rho(1,t|1,t)$ will be a delta function. I guess I need to apply some sort of regularization method to solve the problem.  
Edit2: I have noticed that this indeed Abel type of equation as it becomes unbounded with rate of $(t-t')^{0.5}$ as $t \rightarrow t'$. 
 A: We can write the integral equation as 
$$
\rho(1,t|0, t') = \int_{t'}^{t}dt'' S(t-t'') K(1,t|1,t'') P(t''|t')
$$
where
$$
S(t-t') = \frac{1}{\Sigma(t-t')}
$$
and 
$$
K(1,t|1,t') = \frac{1}{\sqrt{2\pi}}\exp(- \frac{(1-\bar{x}(t|1,t')^2}{2\Sigma^2(t-t')}).
$$
Now, singularity is only in $S$ and $\lim_{t \rightarrow t'} K = 1/\sqrt{2\pi}$. Also we know $\Sigma \rightarrow 0$ simular to $\frac{\tau}{\sigma}(t-t')^{1/2} \rightarrow 0$ as $t \rightarrow t'$. Here, we can use product integration and midpoint method to evaluate approximation of the integral as
$$
\rho(1, t_n|0, t')=\sum_{i=0}^{n-1} w_{ni} K(1,t_n|1,t_{i+\frac{1}{2}})\hat{P}_{t_i+\frac{1}{2}}
$$
where
$$
w_{ni}= \int_{t_i}^{t_{i+1}} S(t_n-u)du.
$$
Thereafter, one can compute the approximation step by step via
$$
\hat{P}_{t_i+\frac{1}{2}}= \frac{\rho(1, t_n|0, t')}{w_{n n-1} K(1,t_n|1,t_{n-\frac{1}{2}})} - \sum_{i=0}^{n-2}  \frac{w_{ni}}{w_{nn-1}}\frac{K(1,t_n|1,t_{i+\frac{1}{2}})}{K(1,t_n|1,t_{n-\frac{1}{2}})} \hat{P}_{t_i-\frac{1}{2}}.
$$

By consulting with [1], we can provide a better approximation of $K(1,t|1,t')P(t|t')$ to improve the result. The next step is the product trapezoidal method which is constructed by approximating $K(1,t|1,t')P(t|t')$ by piecewise linear functions. Thus, we can write our integral equation as
$$
\rho(1, t_n|0, t') = \sum_{j=0}^{n} w_{nj} K(1,t_n|1,t_{j})\hat{P}_{j}
$$
with
$$
w_{n0}= \frac{1}{h} \int_0^h (h - u)S(t_n-u)du 
$$
$$
w_{nj} = \frac{1}{h} ( \int_{t_j}^{t_{j+1}} (t_j - u)S(t_n -u)du + \int_{t_{j-1}}^{t_j} (u- t_{j-1}) S(t_n-u)du)
$$
$$
w_{nn} =  \frac{1}{h} \int_{t_{n-1}}^{t_n} (u- t_{n-1}) S(t_n-u)du
$$
[1] Linz, P. Analytical and Numerical Methods for Volterra Equations. (Society for Industrial and Applied Mathematics, 1987).
