I'm facing a tricky problem. I need to solve a system of 6 differential equations numerically, but I don't have 6 IVP (initial value problem) conditions, instead I have 6 BVP (boundary valye problem) conditions (3 conditions at $x=0$ and 3 conditions at $x=l$, being $l$ the last point to compute, thus $0<x<l$).
So I need to convert the 3 boundary conditions at $x=l$ into 3 at $x=0$ (through the shooting method for example) so I have all 6 IVPs to solve my problem numerically through RK4. So far so good.
I have two problems to solve. The first has the following form:
\begin{cases} \frac{dT}{dx}&=f_1(g) &(1)\\ \frac{dV}{dx}&=f_2(e) &(2)\\ \frac{dM}{dx}&=f_3(V,g) &(3)\\ \frac{dK}{dx}&=f_4(x,M,T) &(4)\\ \frac{dg}{dx}&=f_5(x,T) &(5)\\ \frac{de}{dx}&=f_6(K) &(6)\\ \end{cases}
With the following BVPs:
\begin{cases} T(0)&=c_1 \\ T(l)&=0 \\ V(0)&=c_2 \\ V(l)&=0 \\ M(0)&=c_3 \\ M(l)&=0 \\ \end{cases}
Through the shooting method I need to find $K(0)$, $g(0)$ and $e(0)$. In this case I can trivially find $g(0)$ because (1) and (5) are linearly dependent. To find $K(0)$ and $e(0)$ I shoot, for example:
\begin{cases} K(0)=+1 \text{ and } e(0)=+1 \\ K(0)=+1 \text{ and } e(0)=-1 \\ K(0)=-1 \text{ and } e(0)=-1 \\ K(0)=-1 \text{ and } e(0)=+1 \\ \end{cases}
Then I record the final values for $V(l)$ and $M(l)$ for each case, make an approximating polynomial for each ($f_V=f(K,e)$ and $f_M=f(K,e)$ respectively), and solve the resulting system for $K$ and $e$:
\begin{cases} f_V(K,e)=0 \\ f_M(K,e)=0 \end{cases}
The resulting $K(0)$ and $e(0)$ give me fairly good values (if I use approximating polynomials of order 1 I have an error of $10^-8$, more than acceptable for what I need). Thus, my first problem is solved.
The second problem is the real issue. The system is now:
\begin{cases} \frac{dT}{dx}&=f_1(e,g)\\ \frac{dV}{dx}&=f_2(e,g) \\ \frac{dM}{dx}&=f_3(V,e,g) \\ \frac{dK}{dx}&=f_4(x,M,T) \\ \frac{dg}{dx}&=f_5(x,T) \\ \frac{de}{dx}&=f_6(K) \\ \end{cases}
Also, there are now hyperbolic tangents involved (I have $g(x)$ inside $tanh()$ for example). I can't find any of $K(0)$, $g(0)$ or $e(0)$ trivially because now all equations depend on the others. So, my initial idea was to use the same shooting methodology but creating polynomials with 3 variables ($K$, $g$ and $e$) instead of just $K$ and $e$ like I used before. This isn't producing viable results, even if I use higher orders of polynomials.
My question is: what can I do ? I don't mind changing methodologies or solving using other methods, I just need some directions so I can correctly solve this. I wholeheartedly welcome different perspectives on this. I'm using maxima, for what it's worth.
Also, I can determine the correct $K(0)$, $g(0)$ and $e(0)$ (through an external Maple file), and if I input those into my program I get correct results, so I know all the equations are correctly inputted.
Thank you !
EDIT2:
The system I'm aiming to solve is the following:
\begin{cases} \frac{dT(x)}{dx}&={{3339\,{\it g(x)}\,\tanh \left({{50000\,\sqrt{8427\,{\it g(x)}^2+30036 \,{\it e(x)}^2}}\over{163611}}\right)}\over{\sqrt{8427\,{\it g(x)}^2+ 30036\,{\it e(x)}^2}}}\\ \frac{dV(x)}{dx}&={{12600\,{\it e(x)}\,\tanh \left({{50000\,\sqrt{8427\,{\it g(x)}^2+ 30036\,{\it e(x)}^2}}\over{163611}}\right)}\over{\sqrt{8427\,{\it g(x)}^ 2+30036\,{\it e(x)}^2}}} \\ \frac{dM(x)}{dx}&={{5\,\sqrt{8427\,{\it g(x)}^2+30036\,{\it e(x)}^2}\,{\it V(x)}-10017\, {\it g(x)}\,\tanh \left({{50000\,\sqrt{8427\,{\it g(x)}^2+30036\, {\it e(x)}^2}}\over{163611}}\right)}\over{5\,\sqrt{8427\,{\it g(x)}^2+ 30036\,{\it e(x)}^2}}} \\ \frac{dK(x)}{dx}&=-{{60092431565\,x+15000000000\,{\it T(x)}+25000000000\,{\it M(x)}- 2410745228515}\over{48076923076923}} \\ \frac{dg(x)}{dx}&={{3281046763449\,x+1069000000000\,{\it T(x)}-156626689476919}\over{ 5250000000000000}} \\ \frac{de(x)}{dx}&={\it K(x)} \\ \end{cases}
Boundary conditions are:
\begin{cases} T(0)=200 \\ V(0)=-4.8073945252 \\ M(0)=-47.1403817188 \\ T(l)=0 \\ V(l)=0 \\ M(l)=0 \\ \end{cases}
Where $l=25mm.$
Ok I try to answer but I have no clue if it works. I'll make simpler case where you solve:
$$ u' = f(u,v,x) \\ v'=g(u,v,x)$$
But you only know u(0) and u(l).
Denote grid points $0=x_0 < x_1 < \dots < x_n = l$ with $x_{i+1}-x_i = h$. $u(x_i) = u_i$.
Now using central difference we arrive to system of equations. $$ \frac{1}{2h} \begin{pmatrix} 0 & 1 & 0 & \cdots & 0&0 \\ -1& 0 & 1 & \cdots & 0&0 \\ 0 & -1 & 0 & \cdots & 0&0 \\ \vdots & \vdots &\vdots& \ddots & \vdots &\vdots \\ 0& 0&0& \cdots &0& 1\\ 0&0&0&\cdots &-1 & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\u_3 \\ \vdots \\u_{n-2} \\ u_{n-1} \end{pmatrix} + \frac{1}{2h} \begin{pmatrix} -u(0) \\ 0 \\ 0 \\ \vdots \\0 \\ u(l) \end{pmatrix} =\begin{pmatrix} f(u_1,v_1,x_1) \\ f(u_2,v_2,x_2) \\ f(u_3,v_3,x_3) \\ \vdots \\f(u_{n-2},v_{n-2},x_{n-2}) \\ f(u_{n-1},v_{n-1},x_{n-1}) \end{pmatrix} $$
$$ \frac{1}{2h} \begin{pmatrix} -2 & 2 & 0 & \cdots & 0&0 \\ -1& 0 & 1 & \cdots & 0&0 \\ 0 & -1 & 0 & \cdots & 0&0 \\ \vdots & \vdots &\vdots& \ddots & \vdots &\vdots \\ 0& 0&0& \cdots &0& 1\\ 0&0&0&\cdots &-2 & 2 \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \\v_3 \\ \vdots \\v_{n-1} \\ v_{n} \end{pmatrix} =\begin{pmatrix} g(u_0,v_0,x_0) \\ g(u_1,v_1,x_1) \\ g(u_2,v_2,x_2) \\ \vdots \\g(u_{n-1},v_{n-1},x_{n-1}) \\ g(u_{n},v_{n},x_{n}) \end{pmatrix} $$
When discretizing second equation I used forward and backward difference at point 0 and l respectively.
Now you have system of nonlinear equations which you can try to solve.