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, Γ Here is the code that creates S1: ZeroMatrix = Table[0, {i, 1, 3}, {j, 1, 3}];(*initialize zeros matrix*) x11 = ZeroMatrix; x11[[1, 1]] = 1; x22 = ZeroMatrix; x22[[2, 2]] = 1; x33 = ZeroMatrix; x33[[3, 3]] = 1; x12 = ZeroMatrix; x12[[1, 2]] = 1; x13 = ZeroMatrix; x13[[1, 3]] = 1; x23 = ZeroMatrix; x23[[2, 3]] = 1; x21 = Transpose[x12]; x31 = Transpose[x13]; x32 = Transpose[x23]; Γ1 = Γ; Γ2 = \ Γ;(*defining decoherence rate*) H = Ω (x13 + x31 + x23 + x32) + δ (x22 - x11);(*interaction hamiltonian*) damp32[ρ_] := -I (H.ρ - ρ.H) + Γ1/ 2*(2*(x13.ρ.x31) - (ρ.x33) - (x33.ρ)) + \ Γ2/ 2*(2*(x23.ρ.x32) - (ρ.x33) - (x33.ρ)) + \ Γ3/ 2 *(2 (x11 - x22).ρ.(x11 - x22) - (x11 + x22) ρ - ρ (x11 + x22)); ρ2 = Array[r, {3, 3}];(*Definition of the density matrix*) OL = damp32[ρ2];(*inserting density matrix into master equation*) OLF = Flatten[OL]; ρF = Flatten[ρ2];(*The input and output line vectors*) S1 = Transpose[ Table[D[OLF, ρF[[j]]], {j, 1, 9}]];(*The super operator*) The code plots nothing and also doesn't print any error message. Any suggestions?