clc
clear

L=3;???????????????%?length?of?nozzle
nx=31;?????????????%?number?of?grid?points
dx=L/（nx-1）;???????%?mesh?size??
C=0.5;?????????????%?courant?number
shr=1.4;???????????%?specific?heat?ratio
N=1000;????????????%?number?of?time?steps
x=0:dx:L;??????????%?mesh
A=zeros（1nx）;?????%?section?area
rho=zeros（Nnx）;???%?density
T=zeros（Nnx）;?????%?nondimensional?temperature
V=zeros（Nnx）;?????%?nondimensional?velocity
P=zeros（Nnx）;?????%?nondimensional?pressure
Ma=zeros（Nnx）;????%?Mach?number

%?initial?conditions
for?j=1:nx
????A（j）=1+2.2*（x（j）-1.5）^2;
????rho（1j）=1-0.3146*x（j）;
????T（1j）=1-0.2314*x（j）;
????V（1j）=（0.1+1.09*x（j））*T（1j）^0.5;
????P（1j）=rho（1j）*T（1j）;
????Ma（1j）=V（1j）./（T（1j）.^0.5）;
end


for?i=1:N
????%?time?step
????temp_dt=zeros（1nx）;
????for?ii=1:nx
????????temp_dt（ii）=C*dx/（V（iii）+T（iii）^0.5）;
????end
????dt=min（temp_dt）;???????
????
????%?prediction?step
????d_rho=zeros（2nx-1）;
????d_V=zeros（2nx-1）;
????d_T=zeros（2nx-1）;
????pre_rho=zeros（nx-1）;
????pre_V=zeros（nx-1）;
????pre_T=zeros（nx-1）;
????for?m=1:nx-1
????????d_rho（1m）=-V（im）*（rho（im+1）-rho（im））/dx-rho（im）*（V（im+1）-V（im））/dx-rho（im）*V（im）*（log（A（m+1））-log（A（m）））/dx;
????????d_V（1m）=-V（im）*（V（im+1）-V（im））/dx-（1/shr）*（（T（im+1）-T（im））/dx+（T（im）/rho（im））*（rho（im+1）-rho（im））/dx）;
????????d_T（1m）=-V（im）*（T（im+1）-T（im））/dx-（shr-1）*T（im）*（（V（im+1）-V（im））/dx+V（im）*（log（A（m+1））-log（A（m）））/dx）;
????????pre_rho（m）=rho（im）+d_rho（1m）*dt;????????
????????pre_V（m）=V（im）+d_V（1m）*dt;
????????pre_T（m）=T（im）+d_T（1m）*dt;
????end
????
????%?correction?step
????d_rho_av=zeros（1nx-1）;
????d_V_av=zeros（1nx-1）;
????d_T_av=zeros（1nx-1）;
????for?m=2:nx-1
????????