File contents

clear;clc;

w=0.0;


nr=1000;
r=exp(linspace(-10,log(10),nr));

rho0=0.1*exp(-r.^2);

R0=r;

U0=r;




t=2/3;
rhob=@(x)1./(6*pi*x^2);
rho=rho0+rhob(t);
R=R0;
U=U0;

for i=1:100000
    Rp=gradient(R,r);
    Up=gradient(U,r);
    rhop=gradient(rho,r);
    
    M=cumsum(4*pi*R.^2.*Rp.*rho.*gradient(r));
    
    A=(rho/rhob(t)).^(-w/(1+w));
    
    
    Gam=sqrt(1+U.^2-2*M./R);
    
    Rdot=A.*U;
    rhodot=-A.*rho.*(1+w).*(2*U./R+Up./Rp);
    Udot=-A.*((w/(1+w))*Gam.^2.*rhop./(Rp.*rho)+M./R.^2+4*pi*w*R.*rho);

    dt=0.001;
 
    t=t+dt;
    U=U+Udot*dt;
    R=R+Rdot*dt;
    rho=rho+rhodot*dt;
    subplot(3,1,1)
    loglog(r,rho)
    subplot(3,1,2)
    semilogx(r,U)
    subplot(3,1,3);
    semilogx(r,Gam+U)
    drawnow
end