Hay I am trying to calculate and plot the sensitivity of a solution to some master equation using "ParametricNDSolveValue" but unable to plot the derivative of the solution with respect to the parameter.

y[t_] = {y11[t], y12[t], y13[t], y21[t], y22[t], y23[t], y31[t], 
   y32[t], y33[t]};
equations = 
  Table[(y'[t] - S1.y[t])[[j]] == 0, {j, 1, 
     9}] /. {Ω -> Γ3, Γ -> 
     1};
startConditions = {  y11[0] == 1/3, y12[0] == 0, y13[0] == 0, 
   y21[0] == 0, y22[0] == 1/3, y23[0] == 0, y31[0] == 0, y32[0] == 0, 
   y33[0] == 1/3};
system = Join[equations, startConditions];
solution = 
  ParametricNDSolveValue[system, 
   y33[1], {t, 0, 1}, {δ, Γ3}, 
   MaxSteps -> ∞];
Plot[Evaluate[D[solution[δ, 1], δ]], {δ, -1, 1}]

Where S1 is some $9\times9$ matrix with elements contains the parameters δ, Ω Γ3, Γ

Where S1 is:

S1={{0, 0, I Ω, 0, 0, 0, -I Ω, 0, Γ}, {0, -Γ3 + 2 I δ, I Ω, 0, 0, 0, 0, -I Ω, 0},
{I Ω, I Ω, -Γ + I δ, 0, 0, 0, 0, 0, -I Ω}, {0, 0, 0, -Γ3 - 2 I δ, 0, I Ω, -I Ω, 0, 0},
{0, 0, 0, 0, 0, I Ω, 0, -I Ω, Γ}, {0, 0, 0, I Ω, I Ω, -Γ - I δ, 0, 0, -I Ω},
{-I Ω, 0, 0, -I Ω, 0, 0, -Γ - I δ, 0, I Ω}, {0, -I Ω, 0, 0, -I Ω, 0, 0, -Γ + I δ, I Ω},
{0, 0, -I Ω, 0, 0, -I Ω, I Ω, I Ω, -2 Γ}};

The code plots nothing and also doesn't print any error message. The kernel autoquits when I try to plot.

Whereas:

Plot[Evaluate[solution[δ, 1]], {δ, -1, 1}]

Any suggestions?