How do I input a vector form of ODE's on Runge-Kutta-4 for generality to solve in matlab

Question

I have come across a system of ODE's that are written on vector/Matrix format such that; $Ax'=b$

For simplicity, say the system of ODE's has a vector $x'$ containing first order derivatives of 2-variables: $x$ and $y$ with respect to $t$ such that $x' = [x1' x2']^T $ with given initial conditions on $x1(0)=a0$ and $x2(0)=b0$. Vector $b$ and Matrix $A$ contain dependent variables of $x1, x2, t$ . How would you generally input such a system in Runge-Kutta 4 to solve for $x1$ and $x2$. A general method that will be applicable for a larger ODE system as well. I am using matlab.

For an actual example; suppose

$$\pmatrix{1+x1^2& t-x2\\2x1x2&1}\pmatrix{x1'\\x2'}=\pmatrix{3x1^2x2^2+7x2-4\\x1-x2t+x2}$$

For the given initial conditions $x1(0)=10$ and $x2(0)=20$

How do I input such a format of ODE's into Runge-Kutta-4th order to solve in Matlab, a way that will be general for any larger bigger size of ODEs in a format $Ax'=b$

Thank you very much for your help in advance

Answer 1

If you have procedures A(t,x) and b(t,x) that return the system matrix and rhs vector (this is just as placeholders, one can construct these also in-place), then you can do (using python syntax)

x0 = array((x10,x20))
t = linspace(t0,tf,N+1)

def derivs(t,x):
    dx = linalg.solve(A(t,x),b(t,x))
    return dx

and use the usual RK4 step

def RK4step(f,t,x,h):
    k1 = h*f(t,x)
    k2 = h*f(t+0.5*h,x+0.5*k1)
    k3 = h*f(t+0.5*h,x+0.5*k2)
    k4 = h*f(t+h,x+k3)
    return x+(k1+2*k2+2*k3+k4)/6

etc.

---

Similarly in matlab

x0 = [ x10; x20 ];
t = linspace(t0,tf,N+1);

function dx = derivs(t,x)
    dx = A(t,x)\b(t,x);
end % function

and the usual RK4 step

function xnext = RK4step(f,t,x,h)
    k1 = h*f(t,x);
    k2 = h*f(t+0.5*h,x+0.5*k1);
    k3 = h*f(t+0.5*h,x+0.5*k2);
    k4 = h*f(t+h,x+k3);
    xnext = x+(k1+2*k2+2*k3+k4)/6;
end % function

etc. There is also an ode4 procedure available in the Matlab studios (? Mathworks side with examples and tutorial videos), but it is for scalar examples, needs to be adapted for systems, mainly in growing the results array.