function [delU]=Quasi_1D_exit_bc(delt,delx,Cv,a,To,po,gamma,rhouI,rhouI_1,rhoI,rhoI_1,sI,sI_1)

gamma_fac=(gamma-1)/(gamma+1)   %Used coefficient
uI=rhouI/rhoI;
uI_1=rhouI_1/rhoI_1;

%Get dependent variables
%Temperature
TI=To*(1-gamma_fac*uI^2/a^2);
TI_1=To*(1-gamma_fac*uI_1^2/a^2);
%Pressure
pI=po*(1+gamma_fac*uI^2/a^2)^(gamma/(gamma-1));
pI_1=po*(1+gamma_fac*uI_1^2/a^2)^(gamma/(gamma-1));
%Speed of Sound
cI=sqrt(gamma*pI/rhoI);
cI_1=sqrt(gamma*pI_1/rhoI_1);

%Average the variables
u=(uI+uI_1)/2;          p=(pI+pI_1)/2;
c=(cI+cI_1)/2;          rho=(rhoI+rhoI_1)/2;

%Eigenvalues
lam_1=u*delt/delx;
lam_2=(u+c)*delt/delx;
lam_3=(u-c)*delt/delx;
%Residuals
R1=(-lam_1)/(1+lam_1)*(rhoI-rhoI_1-(1/c^2)*(pI-pI_1));
R2=(-lam_2)/(1+lam_2)*(pI-pI_1+rho*c*(uI-uI_1))-...
   gamma*delt*p*u/(1+lam_2)*((sI-sI-1)/(sI*delx));
R3=(-lam_3)/(1+lam_3)*(pI-pI_1-rho*c*(uI-uI_1))-...
   gamma*delt*p*u/(1+lam_3)*((sI-sI-1)/(sI*delx));

%Use Residuals to obtain differences
if (u/c) < 1
    delp=(R2+R3)/2;
else
    delp=0;
end

delrho=R1+delp/c^2
delu=(R2-delp)/(rho*c)

pI=pI+delp;                 
uI=uI+delu;
rhoI=rhoI+delrho;
TI=pI/((gamma-1)*Cv*rhoI);
eI=rhoI*(Cv*TI+.5*uI^2);

delU=[rhoI;rhoI*uI;eI];