Solving a first-order nonlinear differential equation system

Hi, I wish to solve the following system :

f1( y1(t) , y2(t) , y3(t) ) = 0

f2( y1(t) , y2(t) , y3(t) ) = 0

y3'(t) = f3( y1(t) , y2(t) , y3(t) )

for y1(t), y2(t) and y3(t) where f1, f2 and f3 are arbitrary non linear functions. This is an initial value problem where y1(0), y2(0) and y3(0) are known. In order to use one of the ODE solvers, I have to specify a function of the form dydt=f(t,y), where in my case it would be intuitive to impose y=[y1,y2,y3]. However, I have no explicit formulation for y1'(t) and y2'(t), meaning that I cannot give values for elements dydt(1) and dydt(2) in the derivative function.

I have thought to use a simpler vector of the form y=[y3], with dydt corresponding to y3'(t), and where y1 and y2 would be given as extra parameters. I would use y1, y2 and y3 as initial points to first solve f1 and f2 using fsolve for example, and then use the outcome to compute y3'(t) (i.e. dydt) using the third equations of my system. However, with this approach, I solve a system of 2 equations (f1 and f2) and 3 unknowns (y1,y2,y3). This removes the constraint of the third equation and could lead to poor results. Moreover, iterations are performed automatically by the solver and I cannot see how I would force it to consider as next initial point the solution found at the previous iteration.

Would anyone have a better approach to solve such a problem?