I'm trying to plot phase analysis for the Friedmann's equation, however when running the protocol and trying to plot these phase analysis, it comes up with error that rho wasn't allowed in the set of options!

Here's the code i'm having difficulty with!?!

> SFsys := proc (lambda, W, ics, T, v, u) local sys, flatF, rho0line, plot0, plotDSp, plotDSm, plotE, flatF_Plot, IC, df1, df2;

sys := eval(DSys, {Lambda = lambda, gamma = W});

flatF := eval(solve(subs(K = 0, Feq), rho), Lambda = lambda);

rho0line := plot(0, H = (-1)*v*1.1 .. v*1.1, color = blue, thickness = 3);

flatF_Plot := plot(flatF, H = (-1)*v*1.1 .. v*1.1, 0 .. u, color = blue, thickness = 3);

df1 := dfieldplot(sys, [H(t), rho(t)], t = -T .. T, rho = 0 .. u, H = -v .. v, arrows = MEDIUM,

axes = BOXED, color = grey);

IC := ics; df2 := DEplot(sys, [H(t), rho(t)], t = -T .. T, IC, linecolour = black,

rho0 = 0 .. u, H = -v .. v, stepsize = 0.5e-1, arrows = NONE, thickness = 2,

method = classical[rk4], scene = [H, rho], axes = BOXED);

if lambda = 0 then plot0 := plot([0.2e-1*v*cos(t), 0.1e-1*u*sin(t), t = 0 .. 2*Pi], color = black, thickness = 5); display({rho0line, plot0, flatF_Plot, df1, df2})

elif 0 < lambda then plotDSp := plot([sqrt((1/3)*lambda)+0.2e-1*v*cos(t), 0.1e-1*u*sin(t),

t = 0 .. 2*Pi], color = black, thickness = 5); plotDSm := plot([-sqrt((1/3)*lambda)+0.2e-1*v*cos(t), 0.1e-1*u*sin(t), t = 0 .. 2*Pi], color = black, thickness = 5); if -1/3 < W then plotE := plot([0.2e-1*v*cos(t), 2*lambda/(1+3*W)+0.1e-1*u*sin(t), t = 0 .. 2*Pi], color = black, thickness = 5); display({rho0line, plotDSp, plotDSm, plotE, flatF_Plot, df1, df2}) else display({rho0line, plotDSp, plotDSm, flatF_Plot, df1, df2}) end if elif lambda < 0 and W < -1/3 then plotE := plot([0.2e-1*v*cos(t), 2*lambda/(1+3*W)+0.1e-1*u*sin(t), t = 0 .. 2*Pi], color = black, thickness = 5); display({rho0line, plotE, flatF_Plot, df1, df2}) else display({rho0line, flatF_Plot, df1, df2}) end if end proc;

> ic2 := [[H(0) = VectorCalculus[`-`](1), rho(0) = 1], [H(0) = VectorCalculus[`-`](.5),

rho(0) = 1], [H(0) = 0, rho(0) = 1], [H(0) = .5, rho(0) = 1], [H(0) = 1, rho(0) = 1],

[H(0) = VectorCalculus[`-`](1), rho(0) = VectorCalculus[`-`](1)], [H(0) = VectorCalculus[`-`](.5),

rho(0) = VectorCalculus[`-`](1)], [H(0) = 0, rho(0) = VectorCalculus[`-`](1)], [H(0) = .5,

rho(0) = VectorCalculus[`-`](1)], [H(0) = 1, rho(0) = VectorCalculus[`-`](1)]];

> ic1 := [seq(seq([H(0) = i, rho(0) = j], i = VectorCalculus[`-`](2) .. 2), j = 0 .. 4)];

> ic3 := [seq([H(0) = 0, rho(0) = VectorCalculus[`+`](.2, VectorCalculus[`*`](.5, i))], i = 0 .. 8)];

> ic4 := [seq([H(0) = VectorCalculus[`+`](.2, VectorCalculus[`*`](.5, i)), rho(0) = 12], i = 1 .. 9),

seq([H(0) = VectorCalculus[`+`](VectorCalculus[`-`](.2), VectorCalculus[`*`](.5, i)), rho(0) = 12],

i = VectorCalculus[`-`](9) .. VectorCalculus[`-`](2))];

> ic5 := [seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = VectorCalculus[`-`](9) .. 9)];

> ic6 := [seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = VectorCalculus[`-`](9) .. VectorCalculus[`-`](1)),

seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = 1 .. 9)];

> ic7 := [seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = VectorCalculus[`-`](9) .. VectorCalculus[`-`](5)),

seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = VectorCalculus[`-`](3) .. VectorCalculus[`-`](1)),

seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = 1 .. 3), seq([H(0) = VectorCalculus[`*`](.5, i),

rho(0) = 12], i = 5 .. 9)];

> ic8 := [seq([H(0) = VectorCalculus[`-`](2), rho(0) = VectorCalculus[`*`](2, i)], i = 1 .. 3), seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = VectorCalculus[`-`](3) .. VectorCalculus[`-`](1)), seq([H(0) = VectorCalculus[`*`](.5, i), rho(0) = 12], i = 1 .. 3), seq([H(0) = 2, rho(0) = VectorCalculus[`*`](2, i)], i = 1 .. 3), [H(0) = 0, rho(0) = 6], [H(0) = 0, rho(0) = 9]];

> SFsys(0, 0, ic7, 3, 2, 12);