I am looking to solve (in Python) a differential algebraic equation of the form x'(t) = f(x(t),y(t)) subject to g(x)=0 for a function g:R^n->R^m defining the constraints on the state variable x. I think it should be possible to do that without needing to use fsolve and hence manually building up a numerical method. For example, it seems reasonable this is possible using sci.integrate.solve_ivp or something similar. Can you suggest some specific effective option? Continue reading...
