p = 1e5*ones(80,1);
T = p/100;
M = 29*ones(80,1);
rho = 1.2*ones(80,1); rho0 = rho;
nC = 80;
nF  = 81;
rf = (0:.01:.8)';
r = 0.5*(rf(1:end-1) + rf(2:end));
dr = 0.01;
dt = 1e-7;
uf = zeros(81,1);
rhof0 = interpolate(rho);
uf0 = uf;
URFu = 0.5; URFp = 0.5;

rho = p./T/8314.*M;
rhof = interpolate(rho);

Auf = diag(rhof(2:nF).*rf(2:nF).*(dr/dt + [uf(3:nF) - uf(1:nF-2);uf(nF) - uf(nF-1)]/2));

buf = dr/dt*rf(2:nF).*rhof0(2:nF).*uf0(2:nF) - rf(2:nF).*[p(2:end) - p(1:end-1);0];

%uInitRes = res_compute(Auf,buf,uf(2:end));

uyf(2:nF) = Auf\buf;
uf = (1 - URFu)*uf + URFu*uyf';

mFlux = uf.*rhof.*rf;



ppc = M/8314./T.*r*dr/dt;
ppdiff = rhof.*rf./[1;diag(Auf)];

App = spdiags(ppc + [ppdiff(2:nF-1);0] + ppdiff(1:nF-1),0,nC,nC);

App = spdiags(-ppdiff(1:nF-1),1,App);
App = spdiags(-ppdiff(2:nF-1),-1,App);

bpp = (rho0 - rho)*dr/dt.*r - mFlux(2:end) + mFlux(1:end-1);

pp = App\bpp * URFp;


pInitRes = max(abs(((rho - rho0)/dt*dr.*r + mFlux(2:end) - mFlux(1:end-1))./(max(abs(mFlux)) + 1e-15)));

p += pp;% p + pp;
rho += pp.*M/8314./T;
rhof = interpolate(rho);
uf(2:end-1) -= (pp(2:end)-pp(1:end-1)).*rf(2:end-1)./(diag(Auf)(1:end-1));
uf(end) = uf(end-1)*rf(end-1)/rf(end);

mFlux = uf.*rhof.*rf;
u = 0.5*(uf(1:end-1)+uf(2:end));

