# How to plot coupled Differential equations with changing parameters?

I'm a new user of mathematica. What I'm trying to do is following:

I have two coupled non-linear differential equations. Say

$x'(t)$ = $A− (B +1)x + x^2 y$

$y'(t) = Bx − x^2 y$

What I want is to plot the phase trajectories for 3-5 initial conditions on a single plot and also be able to change the parameter B at will which will show me how the nature of the fixed points change as I change the parameter.

The Problem:

I tried:

sol = NDSolve[{x[b]'[t] ==  1 - (b + 1) x[b][t] + x[t]^2*y[t],
y'[t] == 1*x[t] - x[t]^2*y[t], x[0.1] == 1, y[0.1] == 0}, {x,
y}, {t, 0, 50}]

ParametricPlot[{x[t], y[t]} /. sol, {t, 0.1, 50}]


It works very good and fine. But whenever I' trying to put Manipulate so that I can change initial conditions or it's showing error. Also using ParametricNDSolve is showing some error which I am not figuring out at all. Any help is appreciated. Thank you.

pndsol = ParametricNDSolve[{x'[t] == 1 - (b + 1) x[t] + x[t]^2*y[t],
y'[t] == 1*x[t] - x[t]^2*y[t], x[0.1] == 1, y[0.1] == 0}, {x, y}, {t, 0, 50}, {b}];

ParametricPlot[Evaluate@Table[{x[b][t] , y[b][t]} /. pndsol, {b, 0, 1, .1}], {t, 0, 50},
AspectRatio -> 1] pndsol2 = ParametricNDSolve[{x'[t] == 1 - (b + 1) x[t] + x[t]^2*y[t],
y'[t] == 1*x[t] - x[t]^2*y[t], x[0.1] == d, y[0.1] == e}, {x,
y}, {t, 0, 50}, {b, d, e}];

Manipulate[ParametricPlot[Evaluate@Table[{x[b, d, e][t] , y[b, d, e][t]} /. pndsol2,
{b, 0, 1, .1}], {t, 0, 50}, AspectRatio -> 1, PlotRange -> {{0, 10}, {0, 2}}],
{{d, 1}, 0, 2}, {{e, 0}, 0, 2}] Manipulate[ParametricPlot[Evaluate@({x[b, d, e][t] , y[b, d, e][t]} /. pndsol2),
{t, 0, 50}, AspectRatio -> 1, PlotRange -> {{0, 10}, {0, 2}}, Frame -> True],
{{b, 1}, 0, 1, .1}, {{d, 1}, 0, 2}, {{e, 0}, 0, 2}] I post this for illustration in the event it may be useful. In the following a=1.

Manipulate[
Show[ParametricPlot[f[p, loc[], loc[], t], {t, 0, 50},
PlotRange -> {{0, 6}, {0, 6}}, PlotStyle -> Red,
PerformanceGoal -> "Quality"],
StreamPlot[{1 - (p + 1) x + x^2 y, p x - x^2 y}, {x, 0, 6}, {y, 0,
6}]], {p, 1, 4}, {loc, {0, 0}, {6, 6}, Locator},
Initialization :> (nds =
ParametricNDSolve[{x'[t] == 1 - (b + 1) x[t] + x[t]^2 y[t],
y'[t] == b x[t] - x[t]^2 y[t], x == xi, y == yi}, {x,
y}, {t, 0, 50}, {b, xi, yi}];
f[b_, xi_, yi_, t_] := {x[b, xi, yi][t], y[b, xi, yi][t]} /. nds;)] 