# Difficulties in solving a system of differential equations with a matrix differential equation

I'm trying to solve a system of differential equations that looks like this:

Q = 1;
z[t_] := Sqrt[1/3 (y[t]^2/2 + V[x[t]] + r[t])];
z0 := Sqrt[1/3 (y0^2/2 + x0^(2 n)/(2 n) + r0)];
k = 100*z0;

a[t_] := x[t]^2;
b[t_] := 4*r[t]*(k/g[t]) + y[t];
c[t_] := x[t] + y[t]*g[t];
d[t_] := 1 + 4*r[t];

A[t_] := {{a[t], b[t]},{c[t], d[t]}};
B[t_] := {{0}, {-Sqrt[x[t]^alpha]}};
J[t_] := {{j1[t], j2[t]}, {j3[t], j24[t]}};

sys = {x'[t] == y[t]/z[t],
y'[t] == -3 (1 + Q) y[t] - x[t]/z[t],
r'[t] == -4 r[t] + 3 Q*y[t]^2,
g'[t] == g[t],
J'[t] == -A[t].J[t] - J[t].Transpose[A[t]] + B[t].Transpose[B[t]]};

ICs = {x[0] == 15, y[0] == -0.8, r[0] == 0.005, g[0] == 1,
J[0] == {{0, 1}, {1, 0}}};

soln = NDSolve[{sys, ICs}, {x[t], y[t], r[t], J[t], g[t]}, {t, 0, 10},
Method -> {"StiffnessSwitching"}]


The difference between this and my real code is that, A and J must be 4x4 matrices. I have no idea how to resolve this. Mathematica returns the error:

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0..


But I can't figure out how to write such a problem for 2x2 or 4x4 matrices.

-------- Update ---------

In small versions the suggestion worked great, however for my problem:

    ClearAll["Global*"]
Mp = 1;
gstar = 228.75;
Cr := gstar*(Pi^2/30);
lambda = 1;
V0 := 1;
n = 1;
c = 0.012;

V[x_] := x^(2*n)/(2*n);
dV[x_] := D[V[x], x];
ddV[x_] := D[dV[x], x];
T[t_] := (r[t]/Cr)^(1/4);

Q = 1;
z[t_] := Sqrt[(1/3)*(y[t]^2/2 + V[x[t]] + r[t])];
z0 := Sqrt[(1/3)*(y0^2/2 + x0^(2*n)/(2*n) + r0)];
k = 100*z0;
beta = 1;
alpha = 0;
gamma[t_] := c*x[t]^alpha*T[t]^beta;
gammaT[t_] := c*x[t]^alpha*beta*T[t]^(beta - 1);
gammaPhi[t_] := c*alpha*x[t]^(alpha - 1)*T[t]^beta;

fpsi[t_] := 1 + k^2/(3*1^2*z[t]^2) - y[t]^2/6;
frho[t_] := 1/(6*z[t]^2);
fdphi[t_] := y[t]/6;
fphi[t_] := x[t]/(6*z[t]^2);
gpsi[t_] := (c*x[t]^alpha*T[t]^beta)*z[t]*y[t]^2 -
(k^2/(3*1^2))* ((2*k^2)/(a[t]^2*z[t]^2)) - y[t]^2;

grho[t_] := 4 - (((c*x[t]^alpha*beta*T[t]^(beta -
1))*z[t]*T[t])/(4*r[t]))*y[t]^2 - k^2/(3*1^2*z[t]^2);

gphi[t_] := (-(k^2/(3*a[t]^2*z[t]^2)))*(3*z[t]^2*y[t] + (dV[x]
/.{x -> x[t]})) - z[t]*(c*alpha*x[t]^(alpha -
1)*T[t]^beta)*y[t]^2;

gdphi[t_] := (-(k^2/(3*1^2) + 2*
(c*x[t]^alpha*T[t]^beta)*z[t]))*y[t];

hpsi[t_] := (2*x[t]^3)/z[t]^2 +
((c*x[t]^alpha*T[t]^beta)/z[t])*y[t];

hrho[t_] := ((T[t]*(c*x[t]^alpha*beta*T[t]^(beta -
1)))/(4*r[t]*z[t]))*y[t];

hdphi[t_] := 3 + (c*x[t]^alpha*T[t]^beta)/z[t] + (1/z[t])*((1/2)*
(x[t]^2/2/x[t]^2));

hphi[t_] := k^2/(z[t]^2*1^2) + (ddV[x] /. {x -> x[t]})/z[t]^2 +
((c*alpha*x[t]^(alpha - 1)*T[t]^beta)/z[t])*y[t];


Matrices are formed with these components:

  A[t_] := {{fpsi[t], frho[t], fdphi[t], fphi[t]},
{gpsi[t] + 4*r[t]*fpsi[t], grho[t] + 4*r[t]*frho[t],
gdphi[t] + 4*r[t]*fdphi[t], gphi[t] + 4*r[t]*fphi[t]},
{hpsi[t] + 4*y[t]*fpsi[t], hrho[t] + 4*y[t]*frho[t], hdphi[t] +
4*y[t]*fdphi[t], hphi[t] + 4*y[t]*fphi[t]}, {0, 0, -1, 0}};

B[t_] := {{0}, {(-Sqrt[2*(c*x[t]^alpha*T[t]^beta)*z[t]*
(T[t]/1^3)])*y[t]}, {Sqrt[2*(c*x[t]^alpha*T[t]^beta)*
(T[t]/(a[t]*z[t])^3)]}, {0}};

J[t_] := {{j11[t], j12[t], j13[t], j14[t]}, {j21[t], j22[t],
j23[t], j24[t]}, {j31[t], j32[t], j33[t], j34[t]}, {j41[t],
j42[t], j43[t], j44[t]}};

sys = {x'[t] == y[t]/z[t],
y'[t] == -3 (1 + Q) y[t] - x[t]/z[t],
r'[t] == -4 r[t] + 3 Q*y[t]^2,
a'[t] == a[t],
J'[t] == -A[t].J[t] - J[t].Transpose[A[t]] +
B[t].Transpose[B[t]]};

ICs = {x[0] == 15,
y[0] == -0.8,
r[0] == 0.005,
a[0] == 1,
j11[0] = 0, j12[0] = 0, j13[0] = 0, j14[0] = 0, j21[0] == 0,
j22[0] == 0, j23[0] == 0, j24[0] == 0, j31[0] == 0,
j32[0] == 0, j33[0] == 1 + I*100, j34[0] == -1 + I*100,
j41[0] == 0, j42[0] == 0, j43[0] == 1 + 100, j44[0] == 1};

Block[{n, r0, x0, y0, alpha, beta}, {n, r0, x0, y0, alpha,
beta} = {1, 0.005, 15, -0.8, 0, 1};
soln = NDSolve[{sys, ICs}, Flatten[{x[t], y[t], r[t], J[t],
a[t]}],
{t, 0, 10}, Method -> {"StiffnessSwitching"}]]


 Equation or list of equations expected instead of 0 in the first \
argument.


And I don't know how to solve it. It must be something simple but I can't see it.

-------- UPDATE 2 ----------

The problem was the part where I wrote 'V[x[t]]', instead I directly write x[t]^2/2, and in the parts where dV[x] appears I write x[t] etc. ... Now it's settled. Thanks.

• You seemed to have some misplaced backticks in your code, which I tried to fix. You might want to check my edit.. You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is useful for learning how to format your questions and answers. You may also find the meta Q&A, How to copy code from Mathematica so it looks good on this site, helpful Aug 26, 2023 at 13:40
• The following checks your system and shows that some of the functions are parameters are undefined (at least in the posted code): Solve[sys, D[Flatten@{x[t], y[t], r[t], J[t], g[t]}, t]] /. t -> 0 /. Solve[ICs, Flatten@{x[t], y[t], r[t], J[t], g[t]} /. t -> 0] Aug 26, 2023 at 13:44
• Oh thanks, it was my first question, I'll look at this link, it will be very helpful thanks, Aug 26, 2023 at 13:45
• Yes it's true, when writing the example I forgot to add some things like n = 1; z[t_] := Sqrt[1/3 (y[t]^2/2 + x[t] + r[t])]; x0 = 15; y0 = -0.8; r0 = 0.0005;z0 := Sqrt[1/3 (y0^2/2 + x0^(2 n)/(2 n) + r0)]; But, I believe that this type of system can be solved by writing: NDSolve[sys, D[Flatten@{x[t], y[t], r[t], J[t], g[t]}, t] ] /. t -> 0 /. Solve[ICs, Flatten@{x[t], y[t], r[t], J[t], g[t]} /. t -> 0] Aug 26, 2023 at 13:51
• If those particular value cause a problem for NDSolve, then you should edit the question to include them. Otherwise, I think my answer will show the correct approach. Aug 26, 2023 at 13:53

If I fill in arbitrarily chosen values for the missing parameters and functions, NDSolve computes a solution. Note the use of Flatten, which is not necessary for success, but it produces a more useful output. (Personally, I'd leave off the [t] and solve for {x, y,...}, but I like that form of solution better than the x[t] ->... form.)
Block[{n, r0, x0, y0, alpha, V},