拟一维喷管流动的数值解

%First step
x=zeros（131）;
A=zeros（131）;
p=zeros（131）;
t=zeros（131）;
v=zeros（131）;
n=0;
for s=1:31
x（1s）=n;
A（1s）=1+2.2*（n-1.5）^2;
p（1s）=1-0.3146*n;
t（1s）=1-0.2314*n;
v（1s）=（0.1+1.09*n）*sqrt（t（1s））;
 disp（[‘x‘ num2str（s） ‘=‘ num2str（x（1s））]）
 disp（[‘A‘ num2str（s） ‘=‘ num2str（A（1s））]）
 disp（[‘p‘ num2str（s） ‘=‘ num2str（p（1s））]）
 disp（[‘v‘ num2str（s） ‘=‘ num2str（v（1s））]）
 disp（[‘t‘ num2str（s） ‘=‘ num2str（t（1s））]）
n=n+0.1;
end 
%Second step
DSp=zeros（131）;
DSv=zeros（131）;
DSt=zeros（131）;
trix=0.1;
for i=1:30
DSp（1i）=-v（1i）*（p（1i+1）-p（1i））/trix-p（1i）*（v（1i+1）-v（1i））/trix-p（1i）*v（1i）*（log（A（1i+1））-log（A（1i）））/trix;
DSv（1i）=-v（1i）*（v（1i+1）-v（1i））/trix-1/1.4*（t（1i+1）-t（1i））/trix-t（1i）/（1.4*p（1i））*（p（1i+1）-p（1i））/trix;
DSt（1i）=-v（1i）*（t（1i+1）-t（1i））/trix-（1.4-1）*t（116）*（（v（1i+1）-v（1i））/trix+v（1i）*（log（A（1i+1））-log（A（1i）））/trix）;
 %disp（[‘DSp‘ num2str（i） ‘=‘ num2str（DSp（1i））]）
 %disp（[‘DSv‘ num2str（i） ‘=‘ num2str（DSv（1i））]）
 %disp（[‘DSt‘ num2str（i） ‘=‘ num2str（DSt（1i））]）
end
%Third step
trit=zeros（130）;
c=0.5;
for j=1:30
    trit（1j）=c*trix/（v（1j）+sqrt（t（1j）））;
end
detat=min（trit）;
Lp=zeros（130）;
Lv=zeros（130）;
Lt=zeros（130）;
for k=1:30
    Lp（1k）=p（1k）+DSp（1k）*detat;
    Lv（1k）=v（1k）+DSv（1k）*detat;
    Lt（1k）=t（1k）+DSt（1k）*detat;
    %disp（[‘Lp‘ num2str（k） ‘=‘ num2str（Lp（1k））]）
    %disp（[‘Lv‘ num2str（k） ‘=‘ num2str（Lv（1k））]）
    %disp（[‘Lt‘ num2str（k） ‘=‘ num2str（Lt（1k））]）
end
%Forth step
DSLp=zeros（131）;
DSLv=zeros（131）;
DSLt=zeros（131）;
trix=0.1;
for i=2:30
DSLp（1i）=-Lv（1i）*（Lp（1i）-Lp（1