# problems with NDSolve

I'm having some troubles with Mathematica, I hope you could give me some hints in how to solve them. I'm trying to solve a system of coupled differential equations using NDSolve (the renormalization group equations of the MSSM), everything works fine but when I add some extra terms that should be present also something breaks in the code and the solution is not correct. Because of the structure of the equations, the solution to some of them should not be changed with the presence of this extra terms (the running of adimensional parameters and gauge coup ab[ling constants does not depend on the presence of trilinear soft terms, the ones that introduce the error), some of this terms do not depend on them and the solution should not be the same.

And it does not, if I introduce these new "trilinear" terms at[t] and ab[t], nothing works. I get a solution that is not physically correct. I've checked the physics in depth and it seems that the problem is not there, so I guess there is something wrong with my code. I also looked for typos or similar things and found nothing.

What I did to try to understand what was happening is to solve first some of the equations that do not depend on everything inside the system and then use the solutions I get to solve the other part of the system. When doing that (basically just using the solution from a NDSolve into another NDSolve) I get this error,

"Computed derivatives do not have dimensionality consistent with the initial conditions."

I've seen that these can come when trying to solve a system of differential equations using solutions from the numerical resolution of another system of diff equations (solution of NDSolve into another NDSolve). I don't know how could I do this. Any thoughts?

So basically the problem is the following: I do not know why but when introducing extra equations and contributions to a system of differential equations the code breaks, even the solution of variables that have no dependency at all in this extra contributions behave badly. It seems that the numerical algorithm that mathematica uses breaks.

In order to solve this I try to solve first a subgroup of equations that are only coupled to themselves and whose solution should not be affected by the introduction of these at[t] and ab[t], by doing so I thing I can force the solution of these variables to be correct and then see what happens with the other things. The problem is that I don't know if I'm doing this correctly as I get the error described above.

Pleas if anyone has encountered a similar problem before and knows what could be happening I would be really thankful. Also any tips or tricks for Mathematica are welcome. btw, I use Mathematica 7.

Here is a copy of the code.

Qmax = 35.
tanbeta = 10.
htzero = 173/174 Sqrt[1 + tanbeta^2]/tanbeta
hbzero = 4/174 Sqrt[1 + tanbeta^2]
gaugino = 1.5
sfermions = 1.5
higgsmassparam = 1.5
Exp[Qmax] .174
Xu[t_] := 2 ht[t]^2 ( mhu2[t] + mq2[t] + mu2[t]) + 2 at[t]^2
Xd[t_] := 2 hb[t]^2 (mhd2[t] + mq2[t] + md2[t]) + 2 ab[t]^2
system = NDSolve[{ht'[t] ==
1/(16 \[Pi]^2) ht[
t] (6 ht[t]^2 + hb[t]^2 - 16/3 g3[t]^2 - 3 g2[t]^2 -
13/15 g1[t]^2), ht[0] == htzero,
hb'[t] ==
1/(16 \[Pi]^2) hb[
t] (6 hb[t]^2 + ht[t]^2 - 16/3 g3[t]^2 - 3 g2[t]^2 -
7/15 g1[t]^2), hb[0] == hbzero,
g3'[t] == -3 /(16 \[Pi]^2) g3[t]^3, g3[0] == Sqrt[4 \[Pi] 0.118],
g2'[t] == 1 /(16 \[Pi]^2) g2[t]^3, g2[0] == Sqrt[4 \[Pi]/30],
g1'[t] == 33/5 /(16 \[Pi]^2) g1[t]^3, g1[0] == Sqrt[4 \[Pi]/60],
M3'[t] == -6 /(16 \[Pi]^2) g3[t]^2 M3[t], M3[Qmax] == gaugino,
M2'[t] == 2 /(16 \[Pi]^2) g2[t]^2 M2[t], M2[Qmax] == gaugino,
M1'[t] == 66/5 /(16 \[Pi]^2) g1[t]^2 M1[t],
M1[Qmax] == gaugino}, {g3[t], g2[t], g1[t], ht[t], hb[t], M1[t],
M2[t], M3[t]}, {t, 0, Qmax}]

g3bis[t_] = g3[t] /. system
g2bis[t_] = g2[t] /. system
g1bis[t_] = g1[t] /. system
htbis[t_] = ht[t] /. system
hbbis[t_] = hb[t] /. system
M3bis[t_] = M3[t] /. system
M2bis[t_] = M2[t] /. system
M1bis[t_] = M1[t] /. system

Plot[{g1bis[z], g2bis[z], g3bis[z]}, {z, 0, Qmax}, Frame -> True]
Plot[{M1bis[z], M2bis[z], M3bis[z]}, {z, 0, Qmax}, Frame -> True]


This works fine and the solutions I get are correct, I get the dimensionality problem when running the system below. If I solve everything together, it only works without putting at[t] and ab[t] equations and contributions, when putting them I get the bad behaviour. This bad baheviour translates even to the terms written above that as you can see do not depend on at[t] and ab[t].

systema =
NDSolve[{mhu2'[t] ==
3/16/\[Pi]^2  Xu[t] -
1/16/\[Pi]^2 (
6 g2bis[t]^2 M2bis[t]^2 + 6/5 g1bis[t]^2 M1bis[t]^2),
mhu2[Qmax] == higgsmassparam,
mhd2'[t] ==
3/16/\[Pi]^2  Xd[t] -
1/16/\[Pi]^2 (6 g2bis[t]^2 M2bis[t]^2 +
6/5 g1bis[t]^2 M1bis[t]^2), mhd2[Qmax] == higgsmassparam,
mq2'[t] ==
1/16/\[Pi]^2 (Xu[t] + Xd[t]) -
1/16/\[Pi]^2 (32/3 g3bis[t]^2 M3bis[t]^2 +
6 g2bis[t]^2 M2bis[t]^2 + 2/15 g1bis[t]^2. M1bis[t]^2),
mq2[Qmax] == sfermions,
mu2'[t] ==
2/16/\[Pi]^2 Xu[t] - 1/16/\[Pi]^2 (32/3) g3bis[t]^2 M3bis[t]^2 -
1/16/\[Pi]^2 (32/15) g1bis[t]^2 M1bis[t]^2,
mu2[Qmax] == sfermions,
md2'[t] ==
2/16/\[Pi]^2 Xd[t] - 1/16/\[Pi]^2 (32/3) g3bis[t]^2 M3bis[t]^2 -
1/16/\[Pi]^2 (8/15) g1bis[t]^2 M1bis[t]^2,
md2[Qmax] == sfermions(*,at'[t]==(1/16Pi^2)((at[t](18htbis[t]^2+
hbbis[t]^2-(16/3)g3bis[t]^2-3g2bis[t]^2-(13/15)g1bis[t]^2)+2ab[
t]hbbis[t]htbis[t]+htbis[t]((32/3)g3bis[t]^2M3[t]+6g2bis[t]^2M2[
t]+(26/15)g1bis[t]^2M1[t]))),at[Qmax]==0,ab'[t]==(1/16Pi^2)((ab[
t](18hbbis[t]^2+htbis[t]^2-(16/3)g3bis[t]^2-3g2bis[t]^2-(13/
15)g1bis[t]^2)+2at[t]htbis[t]hbbis[t]+hbbis[t]((32/3)g3bis[t]^2M3[
t]+6g2bis[t]^2M2[t]+(26/15)g1bis[t]^2M1[t]))),ab[Qmax]==
0}, {mhu2[t], mhd2[t], mq2[t], mu2[t], md2[t],at[t],ab[
t]}, {t, 0, Qmax}]

mhu2bis[t_] = mhu2[t] /. systema
mhd2bis[t_] = mhd2[t] /. systema
atbis[t_]=at[t]/.systema
abbis[t_]=ab[t]/.systema
Plot[{mhu2bis[z], mhd2bis[z]}, {z, 0, Qmax}, Frame -> True]


You'd better reduce your problem to draw more attentions (and avoid possible downvotes)… and, there's a reduntant (* in your second sample, and, where's the definition of M1[t], M2[t], M3[t] in the second sample? –  xzczd Aug 2 at 9:41