Question: How to plot phase analysis?

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);

 

Please Wait...