Why is this non-vectorized NDSolve faster than vectorized?

Question

I'm trying to get a feeling for stochastic population models. A simple starting place is a birth-death process whose deterministic limit is the logistic equation. The idea is to model the probability p[n][t] of having n individuals at time t, which results in an infinite set of linear differential equations (truncated at nmax to reduce it to a large but finite set).

I have implemented this two ways: non-vectorized and vectorized. I thought the vectorized implementation would improve the speed, but it is actually much slower. Does anyone know why, or have other tips for speeding this up? The code below is fast enough in both cases, but it gets worse for larger nmax.

(* set up parameters *)
nmax = 120; (* max population size*)
k = 60; (* carrying capacity *)
b0 = 2; (* density-independent birth rate *)
d0 = 1; (* density-independent death rate *)
d1 = (b0 - d0) / k; (* density-dependent death rate *)

Do[
  b[n] = b0*n; (* birth rate *)
  d[n] = (d0 + d1*n)*n; (* death rate *)
, {n, 0, nmax}];

tmax = 10^9;  (* max time *) 

Here's the non-vectorized version:

(* set up equations, ICs, and unknowns *)
eqns = Join[Table[
p[n]'[t] == b[n - 1] p[n - 1][t] - b[n] p[n][t] + d[n + 1] p[n + 1][t] - d[n] p[n][t]
, {n, 1, nmax - 1}],
{p[0]'[t] == d[1] p[1][t],
 p[nmax]'[t] == b[nmax - 1] p[nmax - 1][t] - d[nmax] p[nmax][t]}
];
ics = Table[p[n][0] == If[n == k, 1, 0], {n, 0, nmax}];
unks = Table[p[n], {n, 0, nmax}];

(* solve it *)
AbsoluteTiming[
 sol = NDSolve[Flatten[{eqns, ics}], unks, {t, 0, tmax}][[1]];
]

(* {0.028615, Null} *)

Here's the vectorized version that I thought would be faster:

i[n_] := n + 1; (* add 1 to indices in array *)

(* set up transition matrix *)
a = Normal@SparseArray[Join[
  Table[{i[n], i[n]} -> -b[n] - d[n], {n, 0, nmax}],
  Table[{i[n + 1], i[n]} -> b[n], {n, 0, nmax - 1}],
  Table[{i[n - 1], i[n]} -> d[n], {n, 1, nmax}]
]];

(* set up ICs *)
p0 = Normal@SparseArray[i[k] -> 1, {nmax + 1}];

(* solve it *)
AbsoluteTiming[
  sol2 = NDSolve[{p'[t] == a.p[t], p[0] == p0}, p, {t, 0, tmax}][[1]];
]

(* {0.586452, Null} *)

Both give the same results. For fun, here they are:

LogLogPlot[Evaluate[Table[p[n][t], {n, nmax, 0, -1}] /. sol], {t, 0.001, tmax}]

ListPlot[Table[{n, p[n][10^5]}, {n, 0, nmax}] /. sol, PlotRange -> All]
ListPlot[Table[{n, p[n][10^9]}, {n, 0, nmax}] /. sol, PlotRange -> All]

You can see that the system quickly reaches a quasi-stationary distribution around the carrying capacity of n=60, but over a long timespan reaches the true stationary distribution of everyone dead (p[0]==1).

Update:

I tried out some NDSolve options and found two major improvements for the vectorized code that put it on par with the non-vectorized. First, remove the Normal from the definition of a as per @user21's comment. Then try some different options:

a = SparseArray[
  Join[Table[{i[n], i[n]} -> -b[n] - d[n], {n, 0, nmax}], 
  Table[{i[n + 1], i[n]} -> b[n], {n, 0, nmax - 1}], 
  Table[{i[n - 1], i[n]} -> d[n], {n, 1, nmax}]]];

AbsoluteTiming[
  sol2 = NDSolve[{p'[t] == a.p[t], p[0] == p0}, p, {t, 0, tmax},
    Jacobian -> {Automatic, Sparse -> a["PatternArray"]}
  ][[1]];]

(* {0.043259, Null} *)

AbsoluteTiming[
  sol2 = NDSolve[{p'[t] == a.p[t], p[0] == p0}, p, {t, 0, tmax},
    Jacobian -> {a}
  ][[1]];]

(* {0.16169, Null} *)

AbsoluteTiming[
  sol2 = NDSolve[{p'[t] == a.p[t], p[0] == p0}, p, {t, 0, tmax}, 
    Method -> {"IDA", "ImplicitSolver" -> {"Newton", 
     "LinearSolveMethod" -> {"Band", "BandWidth" -> {1, 1}}}}
  ][[1]];]

(* {0.034921, Null} *)

Looks like the first and third versions are winners. Too bad they don't seem to stack:

AbsoluteTiming[
  sol2 = NDSolve[{p'[t] == a.p[t], p[0] == p0}, p, {t, 0, tmax}, 
  Method -> {"IDA", "ImplicitSolver" -> {"Newton", 
   "LinearSolveMethod" -> {"Band", "BandWidth" -> {1, 1}}}},
  Jacobian -> {Automatic, Sparse -> a["PatternArray"]}
][[1]];]

(* {0.033892, Null} *)