You can do it by remembering how many paths are leading to each node, and when discovering a new node - summarize that number.
For simplicity, let's assume you have regular BFS algorithm, that whenever you use an edge (u,v), calls visit(u,v,k), where:
u - the source node of the edge
v - the target node of the edge
k - the distance from the original source to u
In addition to this, assume you have a mapping d:(vertex,distance)->#paths.
This is basically a map (or 2D matrix) that it's key is a pair of vertex and an integer - distance, and its value is the number of shortest paths leading from the source, to that vertex, with distance k.
It is easy to see that for each vertex v:
d[v,k] = sum { d[u,k-1] | for all edges (u,v) }
d[source,0] = 0
And now, you can easily find the number of shortest paths of length k leading to each node.
Optimization:
You can see that "number of shortest paths of length k" is redundant, you actually need only one value of k for each vertex. This requires some book-keeping, but saves you some space.
Good luck!
The first idea that come to my mind is the next:
Let's name start vertex as s and end vertex as e.
We can store two arrays: D and Q. D[i] is a length of the shortest path from s to i and Q[i] is a number of shortest paths between s and i.
How can we recalculate these arrays?
First of all, lets set D[s] = 0 and Q[s] = 1. Then we can use well-known bfs:
while queue with vertexes is not empty
get v from queue
set v as visited
for all u, there's an edge from v to u
if u is not visited yet
D[u] = D[v] + 1
Q[u] = Q[v]
push u into the queue
else if D[v] + 1 == D[u]
Q[u] += Q[v]
The answer is Q[e].
