Hi

I browsed the forums but could not find a solution to my problem. Ii have a system of differential equations; however, I have the solution to one of the functions as an integral of another function in the system, but since the solutions are all numerical I don't know how to implement this. Here's the system and what I tried at first

#######################

eq1 := diff(W(r), r) = -(1-beta*(W(r)-W0))*(M(r)+4*Pi*r^3*p(r))/(beta*r*(r-2*M(r)));

eq2 := diff(M(r), r) = 4*Pi*r^2*rho(r);

W0 := solve(thetaR = theta0-W0, W0);

ini := evalf({M(r0) = 0, W(r0) = W0});

thetadef := theta(r) = W(r)+thetaR;

pdef := p(r) = subs(thetadef, value((8*sqrt(2)*Pi*(1/3))*(beta/(1-beta*(W(r)-W0)))^(5/2)*ApproximateInt((1+(1/2)*beta*x/(1-beta*(W(r)-W0)))^(3/2)*x^(3/2)*(1-exp(x-W(r)))/(exp(x-theta(r))+1), x = 0 .. W(r), output = sum, partition = 50))):

rhodef := rho(r) = subs(thetadef, value(4*sqrt(2)*Pi*(beta/(1-beta*(W(r)-W0)))^(3/2)*ApproximateInt((1+(1/2)*beta*x/(1-beta*(W(r)-W0)))^(1/2)*(1+beta*x/(1-beta*(W(r)-W0)))^2*x^(1/2)*(1-exp(x-W(r)))/(exp(x-theta(r))+1), x = 0 .. W(r), output = sum, partition = 50))):

vraddef:=vrad(x)=(1+r)^(-1/2)*ApproximateInt((1+r)^(1/2)*M(r)/r^2,r=x..10^(20),output=sum,partition=50):

dsys := evalf(subs(pdef, rhodef, {eq1, eq2, op(ini)})):

dsol:=dsolve(dsys,numeric,output=listprocedure,range=r0..rf,parameters=[theta0,beta,thetaR]);

Wsol:=eval(W(r),dsol):

Msol:=eval(M(r),dsol):

vradsol := subs(M(r) = Msol(r), rhs(vraddef))

vproc1 := proc (X, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(vradsol, [x = X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc

#########################à

this doesn't work. If I try to use the diff equation for vrad:

diff(vrad(r), r) = -(1+r)^(1/2)*M(r)/r^2

ini := evalf({M(r0) = 0, W(r0) = W0, vrad(rf)=0});

dsys := evalf(subs(pdef, rhodef, {eq1, eq2, eq3, op(ini)})):

dsol:=dsolve(dsys,numeric,output=listprocedure,range=r0..rf,parameters=[theta0,beta,thetaR]);

########################

it says it can't solve.

How to deal with this?

Thanks for the help