You mentioned numeric methods, so I'd like to offer a suggestion: First I'd listen to what it's tellin' me: The rate at which the function of two variables is changing in the t direction is the "sum" of some "attenuation" of the function at in the range . The attenuation is of course the composite function . So right off the bat, I'd put this one on the back burner and look at the simpler initial interval problem:

To me, that's conceptually the same and note how I need to know what the function is not at a point as in an ordinary IVP, but along the interval of integration at t=0. Now, if I start at t=0 and keep a running tally of what the function looks like at the previous time period, I can numerically integrate the integral and treat this like an ordinary differential equation and solve it numerically. So, I'd then break up the interval into say for starters, just ten points. Let me pick just one of the points . Now for that point, I solve the ordinary IVP numerically:

via Runge-Kutta. Now, do that ten times at each time period, and always before going to the next time period, numerically integrate a least-square fit of the data points you calculated for the previous time period. This numerical value, then becomes the approximation for the derivative for the next time step.

I'd work on this simple one a while to see what happens and if the results look encouraging, I would then gradually build up the integrand to the form of the original problem. That is, I'd first add the part maybe, say . Then next try a simple form for K such as . Keep playing with it that way until I build up to the problem I need and always be tolerant of the fact that this may not work or the approach may break down at some point.