I have implemented Gray-Scott reaction diffusion in Iron Python 2, mapping grayscale values between 0 and 1. But since it only gets interesting from about 6400 steps, it takes a lot of time.

So can someone smarter / more proficient in Python than me help me make this code more efficient?

Here is the code:

import random as rnd


dA = 1
dB = 0.5
feed = 0.055
kill = 0.062

wHalb = int(width/2)
hHalb = int(height/2)

#decides if gAA is random or 1
if gAA == 0:
    gAA = rnd.random()


def main():
    #create base grid of chemical A and B
    gridA = [ [gAA for i in range(width)] for i in range(height)]
    gridB = [ [0 for i in range(width)] for i in range(height)]
    
    nextA = [ [1 for i in range(width)] for i in range(height)]
    nextB = [ [0 for i in range(width)] for i in range(height)]
    
    color = [ [0 for i in range(width)] for i in range(height)]
    
    
    for x in range (wHalb-baseBlock,wHalb+baseBlock):
        for y in range (hHalb-baseBlock,hHalb+baseBlock):
            gridB[x][y] = 1
    
    x, y, i, j = 0, 0, 0, 0
    
    for n in range(steps):
        for x in range (width):
            for y in range (height):
                a = gridA[x][y]
                b = gridB[x][y]
                nextA[x][y] = (a + (dA * laplaceA(x, y, gridA)) - (a * b*b) + (feed * (1 - a)))
                nextB[x][y] = (b + (dB * laplaceB(x, y, gridB)) + (a * b*b) - ((kill + feed) * b))
        
        tempA = gridA
        gridA = nextA
        nextA = tempA
        tempB = gridB
        gridB = nextB
        nextB = tempB
    color = [ [(nextA[i][j] -  nextB[i][j]) for i in range(width)] for j in range(height)]
    return color


def laplaceA(x, y, gridA):
    sumA = 0;
    xS = x - 1
    xE = x + 1
    yS = y - 1
    yE = y + 1
    
    if (x == 0):
        xS = width-1
    
    if (y == 0):
        yS = height-1
    
    if (x == width - 1):
        xE = 0;
    
    if (y == height - 1):
        yE = 0;
    
    sumA = sumA + gridA[x][y] * -1
    
    sumA = sumA + gridA[xS][y] * 0.2
    sumA = sumA + gridA[xE][y] * 0.2
    sumA = sumA + gridA[x][yE] * 0.2
    sumA = sumA + gridA[x][yS] * 0.2
    
    sumA = sumA + gridA[xS][yS] * 0.05
    sumA = sumA + gridA[xE][yS] * 0.05
    sumA = sumA + gridA[xS][yE] * 0.05
    sumA = sumA + gridA[xE][yE] * 0.05
    return sumA

def laplaceB(x, y, gridB):
    sumB = 0
    
    xS = x - 1
    xE = x + 1
    yS = y - 1
    yE = y + 1
    
    if (x == 0):
        xS = width-1
    if (y == 0):
        yS = height-1
    if (x == width - 1):
        xE = 0
    if (y == height - 1):
        yE = 0
    
    
    sumB = sumB + gridB[x][y] * -1
    
    sumB = sumB + gridB[xS][y] * 0.2
    sumB = sumB + gridB[xE][y] * 0.2
    sumB = sumB + gridB[x][yE] * 0.2
    sumB = sumB + gridB[x][yS] * 0.2
    
    sumB = sumB + gridB[xS][yS] * 0.05
    sumB = sumB + gridB[xE][yS] * 0.05
    sumB = sumB + gridB[xS][yE] * 0.05
    sumB = sumB + gridB[xE][yE] * 0.05
    return sumB

aOut = main()