This is my first question on Stack Exchange - I welcome any suggestions if my approach to asking it does not match the usual conventions around here.
So: I need to solve a system of two second-order, non-linear, coupled differential equations for the functions f and g which depend only on one free variable, i.e. f(x) and g(x). The functions f and g are defined on a finite interval (e.g. [0;10]) and the values of these functions are given at the boundaries (boundary value problem). I need to do this numerically because there is no analytical solution.
I have been scouting a variety of numerical methods and it seems the so-called relaxation method as outlined in ''Numerical Recipes'' by Press et al. might be the way to go. Since I have no numerics background, I may be mistaken in assuming this method is standard knowledge; I will be more than happy to briefly explain what it does if required.
Now, from what I understand, this method works fine for systems of first-order ODEs. In principle, your usual reduction of order does the trick and I can express my system as four coupled first-order ODEs. However, my problem is that I have no boundary conditions for the two additional equations, that is, no BCs for the derivatives of f and g. This is the point where I got stuck.
Questions: 1) Does anyone know if this poses a conceptual problem to the method? How can I circumvent it? 2) Would it be possible to apply the aforementioned relaxation method to the second-order system by 'brute force'?
Thank you very much in advance!
EDIT:
The system of first-order ODEs would look like this:
$$ \frac{\partial f(t)}{\partial t} = u(t) \\ \frac{\partial u(t)}{\partial t} = a(t,f,g) \\ \frac{\partial g(t)}{\partial t} = v(t) \\ \frac{\partial v(t)}{\partial t} = b(t,f,g,v),$$ where $a$ and $b$ are non-linear functions of their arguments (am not showing them here, but of course I know how they look like). $f(0)$, $g(0)$, $f(10)$ and $g(10)$ would be known values. Should anyone know of a method that can solve these for sure - and I don't mind if it is a built-in method of Mathematica, MATLAB or the likes - I would appreciate any suggestion. Also, thanks to everyone who already commented!
There should be no problem in the setup, for a 4-dimensional first-order system you get 4 slots for boundary conditions, which are filled by the values of $f$ and $g$ at both of the boundaries.
What ever you do numerically, you are solving a non-linear system of equations, so convergence to some solution is only guaranteed from initial points close to that solution. Homotopy continuation methods solve that problem by starting each sub-problem at the solution of the previous one, thus close to the new solution. But it can end in fold points, that is, when following the homotopy curve the process would turn back in the opposite direction of the homotopy variable.