Dear all,

I have a question regarding the dsolve procedure in Maple. I'm trying to construct a neutron star model using the Tolman-Oppenheimer-Volkoff equation using a polytropic equation of state (EOS) which requires me to solve the ode system:

diff(rho(r),r)=-(rho(r)*(1+K/(5/3-1)*rho(r)^(5/3-1))+K*rho(r)^(5/3))*(4*Pi*r^3*K*rho(r)^(5/3)+m(r))/(K*5/3*rho(r)^(5/3-1)*r*(r-2*m(r)))

diff(m(r),r)=4*Pi*rho(r)*(1+K/(5/3-1)*rho(r)^(5/3-1)

where I have used the EOS P(r)=K*rho(r)^(5/3) and K is a known constant. rho is the density of the star, m it's mass and P the pressure inside the star.

For the initial conditions I have chosen: rho(10^(-10))=rho_0 and m(10^-10)=0. I have chosen r=10^-10 as the innermost point for the integration, since the differential equation for rho is singular at r=0. rho_0 is the central density of the star.

I solve these equations numerically using:

TOV:=dsolve({ode,ics},numeric,output=listprocedure)

where ode is my system of differential equations and ics are my initial conditions. I need now the radius of the star (R_star), which is the maximum value of r, up until which Maple has carried out the integration.

My problem is, I don't know of any efficient way, to do this. What I'm doing currently is defining a procedure TOVr:=rhs(TOV[1]) and I evaluate it at a very high value of r, for which Maple returns me the error message: "Error, (in TOVr) cannot evaluate the solution further right of ..., probably a singularity". I then use the command TOVr('last') to call the maximum value of r and to store it.

I can use the above method, as long as I'm solving the ODEs only for a few different values of rho_0. But I would like to plot m(R_star) for values of rho_0 ranging from 10^(-14) to 10^(-12) in order to find the value of rho_0, for which I can obtain the maximum value for m(R_star). But this requires me to know the value of R_star for every rho_0 and using the above method is not feasible for say hundred different values of rho_0, since I can't write a loop, because it get's terminated as soon as Maple gives me the first error message.

I was thinking of using perhaps the 'events' command in dsolve, to stop the numeric integration once the value for the pressure drops very low, say below 10^(-46), since the radius at which P(r)=0 defines the stellar surface. I tried using:

TOV:=dsolve({ode,ics},numeric, events=[[K*rho(r)^(5/3)-10^(-46),halt]])

but if I try again to evaluate the solution at a large value of r, I get the above error message, and the integration doesn't get canceled, although the value 10^(-46) is bigger than the value for the pressure I would obtain for R_star using TOVr('last') and Maple shouldn't encounter a singularity.

Am I using the 'events' command wrong? And does somebody know of a more efficient method to obtain the maximum value of a variable after carying out a numerical integration using dsolve?

Sorry for the long post and thank you all.

Hi all,

I'm currently writing a code to solve the equations for a gravitational collapse of a star numerically. I'm using two dimensional Arrays for most of the physical quantities in my finite difference equations, with one index for the spatial dimension (I'm only considering spherically symmetric collapses) and one time index. I wrote two nested loops, where the inner one solves the equations from the centre of the star to the surface (with the exception of one quantitiy, where I have to integrate the other way around) and the outer loop advances the Arrays in time. My problem is, that Maple doesn't seems to evaluate the two dimensional Arrays. Instead for inserting the appropriate value for e.g. A[3,2], Maple just keeps on using A[3,2] which messes up some formulas, since I have to take roots of the numbers. At least that is what seems to be going on when I'm trying to debug the relevant part of the code and I check for the values, with whicih Maple calculates. I'm also using a one dimensional Array to store the timesteps, I'm using in the code and Maple doesn't seem to have problems evaluating them. Is anybody familiar with this problem? Does Maple have trouble with multidimensional Arrays in loops or have I just messed up my code? I have uploaded my code (I hope it worked). I apologize for my code being messy, but it's the first time I'm working with Maple and using numerical methods. The part which troubles me is labeled "# actual code".

Thank you all.

Oppenheimer-Snyder_Maple_3.mw