I'm trying to derive my own Runge Kutta method of order 4 but I'm running into difficulties. Since, it's a 4th order, it needs to satisfy the following requirements: $$ \sum_{i}b_{i} = 1 \;\text{(first order)} \\ \sum_{i}b_{i}c_{i} = \frac{1}{2} \; \text{(second order)} \\ \sum_{i,j}b_{i}a_{i,j}c_{j} = \frac{1}{6}, \sum_{i}b_{i}c_{i}^{2} = \frac{1}{3} \;\text{(third order)} \\ \sum_{i}b_{i}c_{i}^{3} = \frac{1}{4}, \sum_{i,j}b_{i}c_{i}a_{i,j} = \frac{1}{8}, \sum_{i,j}b_{i}a_{i,j}c_{j}^{2} = \frac{1}{12}, \sum_{i,j,k} b_{i}a_{i,j}a_{j,k}c_{k} = \frac{1}{24}\; \text{(fourth order)} $$
So I've chosen $b_{1} = b_{2} = \frac{1}{2}$ and I just need to compute $a_{1,1}, a_{1,2}, a_{2, 1}, a_{2,2}, c_{1}, c_{2}$. To do so, I'm using the following Matlab code:
syms c1 c2 a11 a12 a21 a22;
b1 = (1/2);
b2 = (1/2);
eqn1 = b1*c1 + b2*c2 == (1/2);
eqn2 = b1*(a11*c1 + a12*c2) + b2*(a21*c1 + a22*c2) == (1/6);
eqn3 = b1*(c2^2) + b2*(c2^2) == (1/3);
eqn4 = b1*(c2^3) + b2*(c2^3) == (1/4);
eqn5 = b1*c1*(a11*c1 + a12*c2) + b2*c2*(a21*c1 + a22*c2) == (1/8);
eqn6 = b1*(a11*(c1^2) + a12*(c2^2)) + b2*(a21*(c1^2) + a22*(c2^2)) == (1/12);
eqn7 = b1*(a11*(a11*c1 + a12*c2) + a12*(a21*c1 + a22*c2)) + ...
b2*(a21*(a11*c1 + a12*c2) + a22*(a21*c1 + a22*c2)) == (1/24);
eqn8 = c1 == a11 + a12;
eqn9 = c2 == a21 + a22;
[c1, c2, a11, a12, a21, a22] = solve([eqn1, eqn2, eqn3, eqn4, eqn5, ...
eqn6, eqn7, eqn8, eqn9], [c1, c2, a11, a12, a21, a22]);
However, I keep getting empty solution. Is there an alternative way to do this and compute the remaining terms, either numerically or computationally?
I do not know at what point Matlab fails to solve this system. It is overdetermined, so numerical solution is likely to fail or to converge to an almost-solution with small residual, which nevertheless is far from the correct solution. Symbolic solution is usually implemented with Gröbner bases, this kind of algebra is not the main focus of Matlab.
Manually, one can solve this system by observing its block structure. First you have the quadrature rule that connects the $c_i$ and $b_i$, $$ \sum_{i=1}^s b_ic_i^m=\frac1{m+1} $$ With $s=2$ you can maximally solve a system of $4$ equations (eqns. (0),(1),(3),(4)), which results in the Gauß-quadrature rule of order $2s=4$ with $b_{1,2}=\frac12$, $c_{1,2}=\frac12\pm\frac1{2\sqrt{3}}$.
Next for the determination of the $a_{ij}$ first consider the equations (2),(5),(6),(8),(9) in which they occur linearly. The system (2),(5) \begin{align} \sum_{i=1}^2 b_i\sum_{j=1}^2 a_{ij}c_j&=\frac16\\ \sum_{i=1}^2 b_ic_i\sum_{j=1}^2 a_{ij}c_j&=\frac18\\ \end{align} can be solved for $d_i=\sum_{j=1}^2 a_{ij}c_j$, \begin{align} b_1d_1+b_2d_2&=\frac16\\ b_1c_1d_1+b_2c_2d_2&=\frac18\\\hline (c_2-c_1)b_1d_1&=\frac{c_2}6-\frac18\\ (c_1-c_2)b_2d_2&=\frac{c_1}6-\frac18\\ \end{align} Then with $d_1$ and eqn. (8) the values $a_{1j}$ are fixed, the same for $a_{2j}$ with the equation for $d_2$ and eqn. (9).
It remains to verify the correctness of the eqns. (6) and (7). If carried out in this order of blocks, also Matlab should be able to arrive at the correct method. If everything goes right, this results in exactly one solution, the implicit order $2s=4$ Gauß method. $$ \begin{array}{c|cc} \frac12-\frac{\sqrt3}6 & \frac14 & \frac14-\frac{\sqrt3}6\\ \frac12+\frac{\sqrt3}6 & \frac14+\frac{\sqrt3}6 & \frac14\\ \hline &\frac12&\frac12 \end{array} $$