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 9X9 matrix with elements contains the parameters δ, Ω  Γ3, Γ 

The code plots nothing and also doesn't print any error message.
Any suggestions?