1 | initial version |

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,

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.