function [delU]=Quasi_1D_entrance_SI_bc(delt,delx,Cv,a,To,po,gamma,U1,U2,s1,s2)

gamma_fac=(gamma-1)/(gamma+1);         %Used coefficient
beta=gamma-1;

rho1=U1(1);              rho2=U2(1);
u1=U1(2)/U1(1);          u2=U2(2)/U2(1);
e1=U1(3);                e2=U2(3);

%Get dependent variables
%Pressure
p1=(gamma-1)*rho1*(e1/rho1-.5*u1^2)
p2=(gamma-1)*rho2*(e2/rho2-.5*u2^2);
%p1=po*(1+gamma_fac*u1^2/a^2)^(gamma/(gamma-1))
%p2=po*(1+gamma_fac*u2^2/a^2)^(gamma/(gamma-1));
%Temperature
%T1=To*(1-gamma_fac*u1^2/a^2);
%T2=To*(1-gamma_fac*u2^2/a^2);
T1=p1/(rho1*(gamma-1)*Cv);
T2=p2/(rho2*(gamma-1)*Cv);

%Speed of Sound
c1=sqrt(gamma*p1/rho1);
c2=sqrt(gamma*p2/rho2);

%Average the variables
u=(u1+u2)/2;          p=(p1+p2)/2;
c=(c1+c2)/2;          rho=(rho1+rho2)/2;
alpha=.5*u1^2;

%Eigenvalues
lam_1=u*delt/delx;
lam_2=(u+c)*delt/delx;
lam_3=(u-c)*delt/delx;

drho_du1=(2*u1*po/((gamma-1)*(gamma+1)*Cv*To*a^2))...
       *(1+gamma_fac*u1^2/a^2)^((2-gamma)/(gamma-1));
dp_du1=(2*u1*po*gamma/((gamma+1)*a^2))...
       *(1+gamma_fac*u1^2/a^2)^(1/(gamma-1));
dp_du1=-u1*(gamma-1)*rho1;
M32=-rho*c*(1-lam_3)+delt*gamma*p1*((s2-s1)/(s1*delx));
M33=1-lam_3+delt*gamma*u1*((s2-s1)/(s1*delx));
R3=(-(u-c)/delx)*(p2-p1-rho*c*(u2-u1))-gamma*p1*u1*((s2-s1)/(s1*delx));
%Form the matrix to solve for the change in non-conservative variables
M11=1;                       M12=-drho_du1;                  M13=0;
M21=0;                       M22=-dp_du1;                    M23=1;
M31=0;                       M32=M32;                        M33=M33;

%Invert this matrix to get the change in the non-conservative variables
M=[M11 M12 M13; M21 M22 M23; M31 M32 M33];
R=delt*[0; 0; R3];
delU_p=M^-1*R;

%Multiply by S to get the conservative variables
Si11=1;               Si12=0;                Si13=0;
Si21=u1;              Si22=rho1;             Si23=0;
Si31=alpha;           Si32=rho1*u1;          Si33=1/beta;

Si=[Si11 Si12 Si13; Si21 Si22 Si23; Si31 Si32 Si33];

delU=Si*delU_p;