ClearAll["Global`*"];
ode1 = y'[t] == -Sin[x[t]]/y[t];
ode2 = x'[t] == -Cos[x[t]] (6 Sin[x[t]] Cos[x[t]] + y[t] (b - c (1 + 3*y[t]^2)))/(2*y[t]^3*(b + c (y[t]^2 - 1)));
ode3 = v'[t] == -(b + c*(y[t]^2 - 1))/(4*y[t]*Cos[x[t]]) + 
    Sin[x[t]]/(2*y[t]^2);
    bc = {x[t0] == 0, y[t0] == Br, v[t0] == Log[Dr]};
    Do[tstar = -3 + i/2;
    sols[i] = ParametricNDSolve[{ode1, ode2, ode3, bc}, {x, y, v}, {t, tstar, 0}, {b, c, Br,Dr}];
    data[i] = Table[{Br, Dr} /. FindRoot[{(y[2, c, Br, Dr][0] - 1) /. sols[i], v[2, c, Br, Dr][0] /. sols[i]}, {{Br, 1}, {Dr, 1}}] // Quiet, {c, 0.3, 2.2,.01}];
    lst[i] = Thread[{data[i][[All, 2]], data[i][[All, 1]]}], {i, 0, 4, 1}]
    ListLinePlot[Table[lst[i], {i, 0, 4, 1}], Frame -> True, FrameLabel -> {{"Br", ""}, {"Dr", ""}}, PlotLegends -> Table[Row[{"tstar =", -3. + i/2}], {i, 0, 4, 1}]]

For the above problem, I want to run two more loops for the parameter Br =1 to 3 with any step size and Dr =20 to 40 with any step size and plot them for several points of t0 like above. Also, for 'FindRoot' ,we will now determine parameters b and c rather than Br and Dr. We will have certaintly many points of Br and Dr where solutions will not exist(singularity) and I need to mark those points where I don't have solutions(mark some rectangle or circle or something else or could fill in some boxes if solution exist vs leave them blank if not exist).