Question: How to solve "Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations"?

Hello,

I am trying to plot some movement equations, but the following error keeps happening:

Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations

when I tryed to rewrite the equation the error became:

Error, (in DEtools/convertsys) invalid specification of initial conditions

Can someone help me with those errors?

Here is the code:

restart;

with(LinearAlgebra);
with(linalg);
with(DifferentialGeometry);
with(VariationalCalculus);
with(Tensor);
with(tensor);
with(Tools);
DGsetup([t, x, y, z], M)

g1 := evalDG(N(t)^2*`&t`(dt, dt)-a(t)^2*(`&t`(dx, dx)+`&t`(dy, dy)+`&t`(dz, dz)))

g1inv := InverseMetric(g1)

C1 := Christoffel(g1)

RM1 := CurvatureTensor(C1)

RM1g := CurvatureTensor(g1)

ContractIndices(RM1, [[1, 3]])

ContractIndices(RM1g, [[1, 3]])

RT := RicciTensor(g1, RM1)

ContractIndices(RT, g1inv, [[1, 1], [2, 2]])

S1 := RicciScalar(g1, RM1)

ContractIndices(RT, g1inv, [[1, 1], [2, 2]])

RM1contra := RaiseLowerIndices(g1inv, RM1, [2, 3, 4])

RM1cov := RaiseLowerIndices(g1, RM1, [1])

Kret1 := ContractIndices(RM1contra, RM1cov, [[1, 1], [2, 2], [3, 3], [4, 4]])

eval(simplify(subs(N(t) = 1, Kret1)))

Lag := a(t)^3*N(t)*Kret1

phi = phi(t)

Lam := -(1/2)*a(t)^3*N(t)*((diff(phi, t))^2/N(t)^2+m^2*phi)

Ltot := Lag+Lam

EqELN := EulerLagrange(Ltot, t, N(t))

EqN := eval(simplify(subs(N(t) = 1, EqELN)))

Eqa := subs(N(t) = 1, Ltot)

EqELa := EulerLagrange(EqaD, t, a(t))

tini := 0

with(DEtools)

with(plots)

tfin = 10

(D(a))(tini) = 0

((D@@2)(a))(tini) = 0

((D@@3)(a))(tini) = 0

phi := a^-3

Dphi := 0

(D(N))(tini) = 0

((D@@2)(N))(tini) = 0

((D@@3)(N))(tini) = 0

sys1 := {EqN, a(tini) = 0.1e-2, (D(a))(tini) = 0.1e-2, ((D@@2)(a))(tini) = 0.1e-2, ((D@@3)(a))(tini) = 1}

tfin = 2

p1 := dsolve(sys1, type = numeric, abserr = 1.*10^(-8), relerr = 1.*10^(-8), range = tini .. tfin)

figN := odeplot(p1, [t, a(t)])

sys2 := {EqELa, a(tini) = 0.1e-6, (D(a))(tini) = 0.1e-4, ((D@@2)(a))(tini) = 0.1e-5, ((D@@3)(a))(tini) = 1, ((D@@4)(a))(tini) = 0.1e-5}

p2 := dsolve(sys2, type = numeric, abserr = 1.*10^(-8), relerr = 1.*10^(-8), range = tini .. tfin)

figa := odeplot(p2, [t, a(t)])