NDSolve with piecewise functions defined over endogenously determined intervals

Question

I want to solve a boundary value problem in which differential equations are defined piecewise over endogenous intervals. The unknowns are the trajectories w[t] and mu[t], defined over [tinf,tsup], and some constants (qb,t0,t1).

Differential equations

weq = w'[t] ==
      Piecewise[{  {Qm[t] + mu[t]/G[w[t]] - qb,           tinf <= t <= t0}, 
                   {mu[t]/((1 - G[w[t]]) G[w[t]]),        t0 <= t <= t1}, 
                   {qb + k - Qm[t] + mu[t]/(1 - G[w[t]]), t1 <= t <= tsup}
               }]

mueq= mu'[t] == 
      Piecewise[{  {-g[w[t]] (1-w[t]+1/2 (qb-Qm[t])^2-(mu[t])^2 1/(2(G[w[t]])^2)),        tinf <= t <= t0}, 
                   {-g[w[t]] (1-w[t]+(mu[t])^2 (2 G[w[t]]-1)/(2(G[w[t]])^2(1-G[w[t])^2)), t0 <= t <= t1}, 
                   {-g[w[t]] (1-w[t]+1/2(qb+k-Qm[t])^2-(mu[t])^2 1/(2(1-G[w[t]])^2)),     t1 <= t <= tsup}
               }]

where $t_{inf}$, $t_{sup}$, $G$ and $Q_m$ are known: tinf=3, tsup=4, G[x_]:=1-Exp[-x] and Qm[t_]:=2t-tsup-1

Boundary conditions

transv1=mu[tinf] == (tsup-tinf)G[w[tinf]]
transv2=mu[tsup] == 0;

Other conditions

Given w[t] and mu[t], the constants (qb,t0,t1) are given by the following equations:

eq1 = qb == Qm[t1]-1+mu[t1]/G[w[t1]]
eq2 = qb == Qm[t0]-mu[t0]/(1-G[w[t0]])
eq3 = Integrate[(1-G[w[t]])(Qm[t]-qb)-mu[t],{t,tinf,t0}] == Integrate[G[w[t]](qb+1-Qm[t])-mu[t],{t,t1,tsup}]

Because boundary conditions are at the extreme of the domain ($t_{inf}$ and $t_{sup}$), the solution must be determined by looking at the three intervals $[t_{inf},t_0]$, $[t_0,t_1]$ and $[t_1,t_{sup}]$ simultaneously.

I have tried using ParametricNDSolve to obtain w[t] and mu[t] as functions of (qb,t0,t1), and then solving the system (eq1,eq2,eq3) with FindRoot to find (qb,t0,t1). But with no success.

Any help would be greatly appreciated. Below is my code

Clear["Global`*"];
Qm[t_] := 2 t - tsup - 1;
G[x_] := 1 - Exp[-x];
g[x_] = D[G[x], x];
eq1 = qb == Qm[t1] - 1 + \[Mu][t1]/G[w[t1]];
eq2 = qb == Qm[t0] - \[Mu][t0]/(1 - G[w[t0]]);
eq3 = Integrate[(1 - G[w[t]]) (Qm[t] - qb) - \[Mu][t], {t, tinf, t0}] == Integrate[G[w[t]] (qb + 1 - Qm[t]) - \[Mu][t], {t, t1, tsup}];

transv1 = \[Mu][tinf] == (tsup - tinf) G[w[tinf]];
transv2 = \[Mu][tsup] == 0;

weq = w'[t] == 
  Piecewise[{{Qm[t] + \[Mu][t]/G[w[t]] - qb, 
     tinf <= t <= t0}, {\[Mu][t]/((1 - G[w[t]]) G[w[t]]), 
     t0 <= t <= t1}, {qb + 1 - Qm[t] + \[Mu][t]/(1 - G[w[t]]), 
     t1 <= t <= tsup}}]
\[Mu]eq = \[Mu]'[t] == 
  Piecewise[{{-g[w[t]] (1 - w[t] + 1/2 (qb - Qm[t])^2 - (\[Mu][t])^2/(
        2 (G[w[t]])^2)), 
     tinf <= t <= 
      t0}, {-g[w[t]] (1 - 
        w[t] + (\[Mu][t])^2 (2 G[w[t]] - 1)/(
         2 (G[w[t]])^2 (1 - G[w[t]])^2)), 
     t0 <= t <= 
      t1}, {-g[w[t]] (1 - w[t] + 
        1/2 (qb + 1 - Qm[t])^2 - (\[Mu][t])^2 1/(2 (1 - G[w[t]])^2)), 
     t1 <= t <= tsup}}]

tinf = 3;
tsup = 4;
solw\[Mu] = 
 ParametricNDSolve[{weq, \[Mu]eq, transv1, transv2}, {w, \[Mu]}, {t, tinf, tsup}, {qb, t0, t1}]

(* Out: {w -> ParametricFunction[ <> ], \[Mu] -> ParametricFunction[ <> ]} *)

wint = w /. solw\[Mu][[1]];
\[Mu]int = \[Mu] /. solw\[Mu][[1]];

eq1 = qb == Qm[t1] - 1 + \[Mu]int[qb, t0, t1][t1]/G[wint[qb, t0, t1][t1]];
eq2 = qb == Qm[t0] - \[Mu]int[qb, t0, t1][t0]/(1 - G[wint[qb, t0, t1][t0]]);
eq3 = Integrate[(1 - G[wint[qb, t0, t1][t]]) (Qm[t] - qb) - \[Mu][
      t], {t, tinf, t0}] == Integrate[G[wint[qb, t0, t1][t]] (qb + 1 - Qm[t]) - \[Mu][t], {t, t1, tsup}];

FindRoot[{eq1, eq2, eq3}, {{qb, 2}, {t0, 3.2}, {t1, 3.5}}]

(* Out: ParametricNDSolve::bvdisc: NDSolve is not currently able to solve boundary value problems with discrete variables. And several other error messages*)