@KIRAN SAJJAN modified_symmetry_work.mw in this worksheet includes the command and executes it. you may get expected plots 

colors := [red, blue, green, black]:   # one color per Br value
plots_all1 := []:

for j from 1 to nops(BrVals) do
    
    # Step 1: Update ODE systems for this Br
    odesys[1] := eval(OdeSys[1], {R=0.6, Br=BrVals[j]}):
    odesys[2] := eval(OdeSys[2], {R=0.6, Br=BrVals[j]}):
    
    # Step 2: Find A and B using optimization
    ReturnSolutions := 1:
    opt := Optimization:-NLPSolve(TwoDomainsSolver, -1..1, -1..1):
    AB := op(convert(opt[2][1..2], list)):
    
    # Step 3: Solve both regions
    ReturnSolutions := 2:
    S[1], S[2] := TwoDomainsSolver(AB):

    # Step 4: Create curves for u(y) for both regions, same color
    p_region1 := odeplot(S[1], [y, u[1](y)], y=-1..0, color=colors[j], thickness=2,linestyle=DASH):
    p_region2 := odeplot(S[2], [y, u[2](y)], y=0..1,  color=colors[j], thickness=2):
    
    # Step 5: Combine them into one curve
    plots_all1 := [op(plots_all1), display([p_region1, p_region2])]:
 
od:

# Step 6: Final display
display(plots_all1,
        legend=[seq(cat("Br = ", BrVals[k]), k=1..nops(BrVals))],
        title="Velocity Profiles for Different Br Values",
        labels=["y", "u(y)"],size=[800,700],'axes'= 'boxed'
        );
use this command in the worksheet to get the plots for different values of Br. 

I have included the plot command in the worksheet  it will generate plots for the br values. However, I do not know the command to generate the table values you mentioned in the comment. If any experts are familiar with it, they may be able to provide an answer.