Revision history [back]

Your system of equations has no "easy" analytical (closed-form) solution It might be possible to get one, but I won't delve into it here.

Your question lacks precision on the meaning of your variables and initial numerical values of your functions. Ihave made a number of (likely incorrect) assumptions :

Data:

a,M=function("a,M")
r=var("r")
Eq1=M(r)*a(r)^2*diff(M(r),r)+(M(r)^2*a(r)+6*a(r))*diff(a(r),r)+1/r^2==0
Eq2=a(r)*M(r)*diff(M(r),r)+7*M(r)*r*diff(a(r),r)+2*M(r)*a(r)==0

desolve_system_rk4 is a function integrating systems of differential equations. It has a set of peculiar requirements :

• Equations are expressions of the first-order derivatives of the functions of interest :

.

var(("da,dM"))
Sys1=[u.subs({diff(a(r),r):da, diff(M(r),r):dM}) for u in (Eq1,Eq2)]
Sys2=solve(Sys1,[da,dM])[0]
Sys3=[u.subs({da:diff(a(r),r), dM:diff(M(r),r)}) for u in Sys2]

• It uses representation of functions by variables (deprecated in Sage...)

.

var("va,vM")
Sys4=[u.rhs().subs({a(r):va, M(r):vM}) for u in Sys3]

Now we can solve the system, conditional to our assumptions on the initial values :

NSol=desolve_system_rk4(Sys4,[va,vM],ics=[1,6,0],ivar=r,end_points=exp(3).n())

The result is a list of (independent variable, value1, value2,...) "points". We can plot those points (and functions thereof), use them to create interpolation functions, etc...

line([[log(r),M-r] for r,a,M in NSol])

HTH,