Here is my approach that you may try:

import numpy as np
import matplotlib.pyplot as plt

fig=plt.figure(figsize=(8, 8))
columns = 4
rows = 5
for i in range(1, columns*rows +1):
    img = np.random.randint(10, size=(h,w))
    fig.add_subplot(rows, columns, i)

The resulting image:


(Original answer date: Oct 7 '17 at 4:20)

Edit 1

Since this answer is popular beyond my expectation. And I see that a small change is needed to enable flexibility for the manipulation of the individual plots. So that I offer this new version to the original code. In essence, it provides:-

  1. access to individual axes of subplots
  2. possibility to plot more features on selected axes/subplot

New code:

import numpy as np
import matplotlib.pyplot as plt

w = 10
h = 10
fig = plt.figure(figsize=(9, 13))
columns = 4
rows = 5

# prep (x,y) for extra plotting
xs = np.linspace(0, 2*np.pi, 60)  # from 0 to 2pi
ys = np.abs(np.sin(xs))           # absolute of sine

# ax enables access to manipulate each of subplots
ax = []

for i in range(columns*rows):
    img = np.random.randint(10, size=(h,w))
    # create subplot and append to ax
    ax.append( fig.add_subplot(rows, columns, i+1) )
    ax[-1].set_title("ax:"+str(i))  # set title
    plt.imshow(img, alpha=0.25)

# do extra plots on selected axes/subplots
# note: index starts with 0
ax[2].plot(xs, 3*ys)
ax[19].plot(ys**2, xs)  # finally, render the plot

The resulting plot:

enter image description here

Edit 2

In the previous example, the code provides access to the sub-plots with single index, which is inconvenient when the figure has many rows/columns of sub-plots. Here is an alternative of it. The code below provides access to the sub-plots with [row_index][column_index], which is more suitable for manipulation of array of many sub-plots.

import matplotlib.pyplot as plt
import numpy as np

# settings
h, w = 10, 10        # for raster image
nrows, ncols = 5, 4  # array of sub-plots
figsize = [6, 8]     # figure size, inches

# prep (x,y) for extra plotting on selected sub-plots
xs = np.linspace(0, 2*np.pi, 60)  # from 0 to 2pi
ys = np.abs(np.sin(xs))           # absolute of sine

# create figure (fig), and array of axes (ax)
fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=figsize)

# plot simple raster image on each sub-plot
for i, axi in enumerate(ax.flat):
    # i runs from 0 to (nrows*ncols-1)
    # axi is equivalent with ax[rowid][colid]
    img = np.random.randint(10, size=(h,w))
    axi.imshow(img, alpha=0.25)
    # get indices of row/column
    rowid = i // ncols
    colid = i % ncols
    # write row/col indices as axes' title for identification
    axi.set_title("Row:"+str(rowid)+", Col:"+str(colid))

# one can access the axes by ax[row_id][col_id]
# do additional plotting on ax[row_id][col_id] of your choice
ax[0][2].plot(xs, 3*ys, color='red', linewidth=3)
ax[4][3].plot(ys**2, xs, color='green', linewidth=3)


The resulting plot: