Hello, 

I am trying to plot 2 equations divided in 3 intervals (tini..tr, tr..te, te..tfin1), but when i try to solve the first equation it appears:

"Warning, cannot evaluate the solution further left of .66099564e18, probably a singularity"

Is it possible avoid this point?

About the second equation, when i try to solve between "te and tr" the program says that the argument is invalid.

Can someone help me?

Here is the program:

restart;

Digits := 20;

M := 0.122e20*0.152e25;
kappa := evalf(sqrt(8*Pi/M^2));
rrhoCC := 0.769e-41*(0.1e-28*0.561e24)*0.152e25^4;

rrhoM0 := 0.1e51;
rrhoR0 := 0;
rrhoM := proc (t) options operator, arrow; rrhoM0/a(t)^3 end proc;
rrhoR := proc (t) options operator, arrow; rrhoR0/a(t)^4 end proc;
pM := proc (t) options operator, arrow; 0 end proc;
eq1 := diff(a(t), t) = a(t)*kappa*sqrt(rrhoM(t)+rrhoCC)/sqrt(3);
eq2 := 2*(diff(a(t), t, t)) = -(1/3)*kappa^2*(rrhoM(t)+3*pM(t)-2*rrhoCC)*a(t);
tini := 0.1e19;
sys1 := {eq1, a(tini) = 1.0};
sys2 := {eq2, a(tini) = 1.0, (D(a))(tini) = 0.284e-17};
with(DEtools);
with(plots);
tfin1 := 0;
te := 10^(-32);
tr := 10^12;
singularities(eq1, a(t));
singularities(eq2, a(t));
p1a := dsolve(sys1, type = numeric, abserr = 10^(-15), relerr = 10^(-15), range = tini .. tr);
p1b := dsolve(sys1, type = numeric, abserr = 10^(-10), relerr = 10^(-10), range = tr .. te);
p1c := dsolve(sys1, type = numeric, abserr = 10^(-10), relerr = 10^(-10), range = te .. tfin1);
fig1a := odeplot(p1a, [t, a(t)]);
fig1b := odeplot(p1b, [t, a(t)]);
fig1c := odeplot(p1c, [t, a(t)], color = red);
display(fig1a, fig1b, fig1c);
p2a := dsolve(sys2, type = numeric, abserr = 1.*10^(-12), relerr = 1.*10^(-12), range = tfin1 .. te);
p2b := dsolve(sys2, type = numeric, abserr = 1.*10^(-12), relerr = 1.*10^(-12), te .. tr);
p2c := dsolve(sys2, type = numeric, abserr = 1.*10^(-12), relerr = 1.*10^(-12), range = tr .. tini);
fig2a := odeplot(p2a, [t, a(t)], color = green);
fig2b := odeplot(p2b, [t, a(t)], color = green);
fig2c := odeplot(p2c, [t, a(t)], color = green);
display(fig2a, fig2b, fig2c);
display(fig1a,fig2a,fig1b,fig2b,fig1c,fig2c);

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)])