Split Boundary Value Problems with Algebraic Equations

Question

Is it possible at all to solve with NDSolve (or other built in function) a split boundary value problem with algebraic equations?

Please look at the following example:

s = NDSolve[
   {
    x'[t] == y[t] + x[t] y[t] + z[t],
    z'[t] == 2*y[t],
    2 x[t] + y[t] - z[t] == 1,
    x[0] == 0,
    z[0] == 0}
   , {x, y, z}, {t, 0, 1}];

This example solves well, and outputs Interpolation Functions.

I compute the final value of variable z:

In[50]:= z[1] /. s
Out[51]:= {0.694658}

When I try the same example with split conditions:

s1 = NDSolve[
   {
    x'[t] == y[t] + x[t] y[t] + z[t],
    z'[t] == 2*y[t],
    2 x[t] + y[t] - z[t] == 1,
    x[0] == 0,
    z[1] == 0.694658}
   , {x, y, z}, {t, 0, 1}];

I get the following error:

NDSolve::bvdae: Differential-algebraic equations must be given as initial value problems.

I am trying to solve a complex problem with this structure, is there any way to tackle it in mathematica?

Edit: Unlike this Minimum Working Example, in the original problem I am trying to solve the variable y cannot be eliminated from the system

Answer 1

Indeed, NDSolve cannot solve this equation as written. However, it is easy enough to eliminate y from the system.

{x'[t] == y[t] + x[t] y[t] + z[t], z'[t] == 2*y[t]} /. y[t] -> 1 + z[t] - 2 x[t]

and then solve and plot

s1 = NDSolve[{Derivative[1][x][t] == 1 - 2 x[t] + 2 z[t] + x[t] (1 - 2 x[t] + z[t]), 
    Derivative[1][z][t] == 2 (1 - 2 x[t] + z[t]), x[0] == 0, 
    z[1] == 0.694658}, {x, z}, {t, 0, 1}]    
Plot[Evaluate[{x[t], 1 + z[t] - 2 x[t], z[t]} /. s1], {t, 0, 1}, 
    AxesLabel -> {t, "x, y, z"}]

Addendum

A similar problem is 63824. Its approach can be applied as follows.

s3 = ParametricNDSolve[{x'[t] == y[t] + x[t] y[t] + z[t], z'[t] == 2*y[t], 
    2 x[t] + y[t] - z[t] == 1, x[0] == 0, z[0] == a}, {x, y, z}, {t, 0, 1}, {a}];
f[w_?NumericQ] := z[w][1] /. s3
w0 = w /. FindRoot[f[w] == 0.694658, {w, .1}]
(* -2.47361*10^-6 *)
Plot[Evaluate[{x[w0][t], y[w0][t], z[w0][t]} /. s3], {t, 0, 1}, 
    AxesLabel -> {t, "x, y, z"}]

which gives the same plot shown above. Note that the value of w0 is zero to roundoff.