%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