I have a large system of non-linear differential equations of the form:
$x_{i}''(y) = f_{i}(x_{1}, x_{2},\ldots,x_{n},x_{1}', x_{2}',\ldots,x_{n}', x_{1}'', x_{2}'',\ldots,x_{i-1}'',x_{i+1}'',\ldots,x_{n}''), \quad i=1,\ldots,n$
which I need to solve numerically. However, the presence of the second-order derivatives in the right-hand side poses a problem. I know I can try and work out how to convert the system to a proper form, but given the amount and complexity of the actual equations, it might prove extremely hard/impossible.
What alternatives do I have? I need only a numerical solution, so I was thinking of perhaps using Newton's method to solve the non-linear system and find the values of $x_{i}''(y)$, but perhaps there are other methods for such problems which I don't know of?
Your system is at least an implicit ODE and might be unsolvable in the ODE framework. The more general framework is that of differential-algebraic equations.
So if there is nothing like the analog of diagonal dominance to save the day, you will need the partial derivatives and Newton method to evaluate the vector of second derivatives, if possible.
If not, then you need to compute the index of the DAE and a reduction to at least an index-1-system.