How should terms be scaled by finite dx and dt in numerical integration of 1D diffusion? I am familiar with numerically integrating systems of ordinary-differential equations, but I feel that I am missing something important in terms of how numerically integrating ODEs differs from numerically integrating PDEs. 
Specifically, I am confused by the fact that in my python program (attached below), when I increase the spatial resolution beyond a certain point, the integration "explodes". 
When integrating ODEs, increasing the temporal resolution is always a better approximation of the "true" system, as a smaller time step (dt) is a better approximation for an infinitesimal dt. But the same does not appear to be the case in my program when I increase the spatial resolution. In this case, too low a resolution is not good, but oddly, too high a resolution is also not good.
Specifically, when dx is much less than 1.0, the integration explodes.  I think that this is either because (1) I am using dx incorrectly in my integrator, or (2) there is something more fundamental that I am missing in my transition from integrating ODEs to integrating PDEs.
The code for calculating the derivative is the following :
next_y[1:-1] += dt*k_r*(y[:-2] - 2.0*y[1:-1] + y[2:])/(dx**2)

it is equivalent to $$y^{p}_{t+1} = y^p + dt \cdot k (-2y^{p}+y^{p-1}+y^{p+1}) / dx^2,$$ where the superscript value refers to discrete spatial index. 
Put into words, this equation takes the difference between each node and each of its neighbours, and scales it by k, the diffusion rate constant, dt the time-step and $dx^2$. I divide by $dx^2$ because the portion inside the parentheses approximates the spatial derivative twice. I am fairly confident that it is this use of dx that is incorrect. But I do not know what the correct usage is.
I have also considered / experimented with scaling the equation by $dx^2$ like this (multiplying by $dx^2$ instead of its reciprocal):
next_y[1:-1] += dt*k_r*(y[:-2] - 2.0*y[1:-1] + y[2:])*(dx**2)

But this makes the system behave qualitatively differently for slightly different values of dx, even when they are both very small.

Question: What is the correct use of dx in this equation? Ideally, I'd love to understand conceptually, and not just be given the right answer!

The code below performs two integrations of the same system, animating the result. 
It should be pretty self-explanatory. The diffusion() function performs the integration, and the animateTwoHistories() function presents an animation of the system changing over time, allowing for an easy comparison of the implications of changes in parameters.
from pylab import *
import matplotlib.animation as animation

def diffusion(# time
        duration = 10.0,
        dt = 0.01,
        #space
        width = 100.0, # cm
        dx = 1.0,
        k_r = 50.0) :

    ## number of discretizations in space and time
    n_x = int(width/dx)
    n_its = int(duration/dt)

    ## initial conditions, all zero but a block of 1s in the middle
    y = zeros(n_x,'f')    
    block_sx = 0.4*width
    block_ex = 0.6*width
    y[int(block_sx/dx):int((block_ex)/dx)] = 1.0

    t = []
    x = linspace(0,width,n_x)    
    h_y = [] ## history of all y values, row=iteration, col=x (spatial location)    
    for it in xrange(n_its) :
        t.append(it*dt)
        h_y.append(y)
        next_y = array(y) # duplicate previous state
        next_y[1:-1] += dt*k_r*(y[:-2] - 2.0*y[1:-1] + y[2:])/(dx**2)
        y = next_y 

    return t,x,h_y

def animateTwoHistories(t1,x1,y1, # data of history one
                        t2,x2,y2  # data of history two
) :
    def update_line(num, data1, data2, l0, l1):
        l0.set_ydata(data1[min(num,len(data1))])
        l1.set_ydata(data2[min(num,len(data2))])
        return l0,l1

    fig1 = plt.figure()   
    l0, = plt.plot(x1, x1, 'k-x')
    l1, = plt.plot(x2, x2, 'r-x')
    ylim(0,1)

    numFrames = max(shape(y1)[0],shape(y2)[0])
    line_ani = animation.FuncAnimation(fig1, update_line, numFrames,
                                       fargs=(y1,y2, l0,l1),
                                       interval=1, blit=False)
    plt.show()

#### changing spatial resolution
t1,x1,y1 = diffusion(dx=0.8) # black plot (fails)
t2,x2,y2 = diffusion(dx=2.0) # red plot (works)

animateTwoHistories(t1,x1,y1,t2,x2,y2)

 A: Looks like it fails when dt*k_r/(dx**2) >= 0.5 (and almost works at the boundary - try k_r = 32)
The discrete model is somewhat like an IIR filter (current state depending linearly on previous state, recursively) and I suspect that is the point where the filter becomes unstable.
Put it another way: find the state transition for y, as an n_x x n_x matrix (it will have zeros except near the diagonal). Then find its eigen decomposition. If any of the eigenvalues have a magnitude >=1, the model will be unstable (because raising it to power $n$ will give a matrix with the same eigen decomposition, but with all the eigenvalues raised to the power $n$; and stability means that it should tend to zero when raised to a high enough power).
[In the code, the two endpoints are never modified; so the 'transition matrix' should actually be the n_x-2 by n_x-2 matrix which updates all the other elements from their previous values, assuming the endpoints are zero. Otherwise the decomposition will always show two eigenvalues which are always 1].
It may be that with this kind of model you want to be quite well into the 'stable' zone in order to ensure that the discretization doesn't compromise the accuracy (not my area of expertise, but that feels right ... ). And as you've shown, that can depend on the model parameters. The eigen decomposition is a good general way to test for that, probably in this case there's a way to solve for stability more directly (and thus confirm or refute my experimental result above). Once known, you can balance dx and dt as needed to make it work for a given  k_r.
A: I think that my problem may have simply been the too high value of the diffusion rate constant k_r. Reducing this from 50.0 to 2.0 seems to make the program work as desired.
