Solving a system of differential equations

This is not a solution. Just a slightly elaborate approach.

Let us consider, for the explanation of the idea, an (over)simplification :

• all constants are $0$, and as well, (see our discussion)

• $\partial J/\partial x=0$ (considered as a constant).

If $A:=\pmatrix{3&4&&&\\1&4&3&&\\&4&3&1&\\&&3&4&5\\3&&4&&1}$, the solution is $\pmatrix{F(t)\\G(t)\\H(t)\\I(t)\\J(t)}=exp(tA).\pmatrix{F(0)\\G(0)\\H(0)\\I(0)\\J(0)}$

In Matlab (I rely on my memory), it would be implemented more or less like that (using "expm" function, not "exp"):

A=[3 4 ...; ... ; 3 0 4 0 1];

V0=[1,3,0,0,0]'; % or any other set of initial values.

n=1000;

T=zeros(5,n);

for k=1:n

t=k/100; % time

V=expm(t*A)*V0;% V0 has to be a column vector

T(1:5,k)=V;

end

plot((1:n)/100,T(1,:)) ; % example : the evolution of variable F

Important remark: The eigenvalues of $A$ are

$$\{8.2, 3.69 \pm 0.8i, -0.29 \pm 1.48i\}$$

It attracts attention on the fact that there is a trend to "global" exponential growth (the dominant eigenvalue $8.2$ is $>0$) with periodic fluctuations.

Edit

One can improve the computational efficiency without impairing the overall precision by using the fact that :

$$\tag{1}\forall t1,t2, \ \text{expm}(t1*A)*\text{expm}(t2*A)=\text{expm}((t1+t2)A)$$

(Remark: it is not true in general that expm$(M_1)*\text{expm}(M_2)=\text{expm}(M_1+M_2)$; but it is true in particular when $M_1$ and $M_2$ can be expressed as polynomials in the same matrix, as is the case here).

If relationship (1) is applied with $t_2=t$ and $t_1=\Delta t$, this gives rise to the fact that it suffices to left-multiply the current state $\text{expm}(t*A)*V_0$ by a multiplicative factor, always the same : $M=\text{expm}((\Delta t)*A)$ to obtain the new current state.

This gives rise to the following modification of the heart of the previous program:

V = V0;

Deltat = 0.01;

M=expm(Deltat*A);

for k=1:n

T(1:5,k)=V;

V=M*V;

end