Solving a "large" differential equation system

edit

Hello,

As said in the title, i'm having issues solving a system involving 7 differential equations.

The variables and initial conditions are defined as follow:

var('g1 g2 g3 lambd y_t y_b y_tau mu g_10 g_20 g_30 lambda_0 y_t0 y_b0 y_tau_0 mu_top')

#constants and ics:

g_10 = (5/3)**(1/2)*0.35830
g_20 = 0.64779
g_30 = sqrt(4*pi*0.1184)
lambda_0 = 0.12604
y_t0 = 0.94018
y_b0 = sqrt(2)*(4/246)
y_tau_0 = sqrt(2)*(1777/246)
mu_top = 173.34

k = (4*pi)**(-1)

Then the functions:

g1 = function('g_1')(mu)
g2 = function('g_2')(mu)
g3 = function('g_3')(mu)

lambd = function('lambda')(mu)

yt = function('y_t')(mu)
yb = function('y_b')(mu)
ytau = function('y_tau')(mu)

And the equations:

dg1 = diff(g1**2,mu,1) == k**2 * g1**4 * 4.1
dg2 = diff(g2**2,mu,1) == k**2 * g2**4 * (-19/6)
dg3 = diff(g3**2,mu,1) == k**2 * g3**4 * (-7)

dlambda = derivative(lambd,mu,1) == k**2 * (lambd*(12*lambd + 6*yt**2 + 6*yb**2 + 2*ytau**2 - (9/2)*g2**2 - (9/10)*g1**2) -3*(yt**4 + yb**4 + ytau**4) + (9/16)*g2**4 + (27/100)*g1**4 + (9/40)*g2*g2*g1*g1)

dyt = diff(yt**2,mu,1) == k*k*yt*yt * ((9/2)*yt**2 + (3/2)*yb**2 + ytau**2 - 8*g3**2 - (9/4)*g2**2 - (17/20)*g1**2) 
dyb = diff(yb**2,mu,1) == k*k*yb*yb * ((3/2)*yt**2 + (9/2)*yb**2 + ytau**2 - 8*g3**2 - (9/4)*g2**2 - (1/4)*g1**2);
dytau = diff(ytau**2,mu,1) == k*k*ytau*ytau * (3*yt**2 + 3*yb**2 + (5/2)*ytau**2 - (9/4)*g2**2 - (9/4)*g1**2)

I tried the following solving methods (of course without the # before the function I used):

des1 = [dg1, dg2, dg3, dlambda, dyt, dyb, dytau]; show(des1)
ics1 = [mu_top, g_10, g_20, g_30, lambda_0, y_t0, y_b0, y_tau_0] ; show(ics1)
dvars1 = [g1,g2,g3,lambd,yt,yb,ytau]; show(dvars1)



#sol1 = desolve_odeint([dg1, dg2, dg3, dlambda, dyt, dyb, dytau], [mu_top, g_10, g_20, g_30, lambda_0, y_t0, y_b0, y_tau_0], times=srange(10,100,0.2), dvars=[g1,g2,g3,lambd,yt,yb,ytau],ivar=mu);show(sol1)

#sol1 = desolve_odeint(des1, ics1, times=srange(10.,100.,1.), dvars=dvars1, ivar=mu)
#sol1 = desolve_system_rk4(des1, vars=dvars1,ics=ics1,ivar=mu,end_points=[0,100]) ;show(sol1)
#sol1 = desolve_system(des1, vars=dvars1,ics=ics1,ivar=mu) ;show(sol1)
#sol1 = desolve_tides_mpfr(des1, ics=ics1, initial=10.3, final=100, delta=0.2, tolrel=10**(-16), tolabs=10**(-16), digits=50)

Unfortunately, none of them seemed to work properly, and i can't find what i did wrong. It'd be great if one of you could help me find my mistakes.

Thanks in advance