function?myjob（）
clcclear
tic
indata
z=15;
n=8;
h=0.01;
Y=[0?20?z?0?z?0?0?0];??%初始值
fid=fopen（‘result.txt‘‘w‘）;
%將數(shù)據(jù)寫入文本文件
for?i=1:n
????fprintf（fid‘%8.4f?‘Y（i））;
end
old4=YFZ（Y）;
Y=RK（nhY）;
%write
for?i=1:n
????fprintf（fid‘%8.4f?‘Y（i））;
end
old3=YFZ（Y）;
Y=RK（nhY）;
%write
for?i=1:n
????fprintf（fid‘%8.4f?‘Y（i））;
end
old2=YFZ（Y）;
Y=RK（nhY）;
%write
for?i=1:n
????fprintf（fid‘%8.4f?‘Y（i））;
end
old1=YFZ（Y）;
while?Y（7）>=0
????[Yold1old2old3old4]=ADM（Yhold1old2old3old4）;
????%write
????for?i=1:n
????????fprintf（fid‘%8.4f?‘Y（i））;
????end
end
fclose（fid）;
Y
Y（1）
%結(jié)果處理
load?result.txt
s=size（result2）;
da=reshape（result8s/8）;
da=da‘;
plot（da（:6）da（:7））title（‘彈道曲線‘）xlabel（‘x/m‘）ylabel（‘y/m‘）
figure
plot（da（:1）da（:2））title（‘速度曲線‘）xlabel（‘t/s‘）ylabel（‘V/（m/s）‘）
a=da（:5）-da（:3）;
figure
plot（da（:1）a）title（‘攻角曲線‘）xlabel（‘t/s‘）ylabel（‘a(chǎn)/rad‘）
toc

%阿當姆斯預估校正法
function?[Yold1old2old3old4]=ADM（Yhold1old2old3old4）
y=Y+h*（55*old1-59*old2+37*old3-9*old4）/24;
f=YFZ（y）;
Y=Y+h*（9*f+19*old1-5*old2+old3）/24;
old4=old3;
old3=old2;
old2=old1;
old1=YFZ（Y）;

%四階龍格庫塔法
function?Y=RK（nhY）
b=[h/2?h/2?h?h];
dy=YFZ（Y）;
oldy=Y;
Y1=[];
for?j=1:3
????for?i=1:n
????????Y1（i）=oldy（i）+b（j）*dy（i）;
????????Y（i）=Y（i）+b（j+1）*dy（i）/3;
????end
????dy=YFZ（Y1）;
end
Y=Y+b（1）*dy/3;

%右端項函數(shù)
function?dy=YFZ（Y）
ts=3.6;
tf=ts+0.5;
T=273.5+20;
Ma=Y（2）/（20*sqrt（T））;
m0=819.06;
S=0.308;
D=0.36;
RHO=1.225;
g=9.8;
q=RHO*Y（2）*Y（2）/2;
if?Y（1）????mro=0;
else
????mro=152;
end
if?Y（1）????P=76930;
????mc=28.224;
elseif?Y（1）>=ts?&&?Y（1）????P=0;
????mc=0;
elseif?Y（1）>=tf?&&?Y（1）<145.15
????P=2270;
????mc=1.094;
else?
????P=0;
????mc=0;
end
m=m0-mro-Y（8）;
a=Y（5）-Y（3）;
%插值子程序的調(diào)用:cx、cya、mza、mzw、Jz
[cxcyamzamzwJz]=inter（Y（1）Maatf）;
dy（1）=1;
dy（2）=P*cos（a/57.3）/m-cx*q*S/m-g*sin（Y（3）/57.3）;
dy（3）=57.3*（P*sin（a/57.3）+cya*a*q*S-m*g*cos（Y（3）/57.3））/（m*Y（2））;
dy（4）=57.3*（mza*a+mzw*Y（4）*D/（57.3*2*Y（2）））*q*S*D/Jz;
dy（5）=Y（4）;
dy（6）=Y（2）*cos（Y（3）/57.3）;
dy（7）=Y（2）*sin（Y（3）/57.3）;
dy（8）=mc;

%插值子程序
function?[cxcyamzamzwJz]=inter（tMaatf）
load?data;
A=[0?1?2?3?4?5?6];
MC=[0.0?0.1?0.3?0.5?0.7?0.8?0.9?0.95];
T=[4.1?28.08?56.29?90.14?145.15];
MM=[0.7?0.75?0.8?0.85?0.9?0.95?1.00];
MCA=[0?1?2?3?4?6?8];
if?t????Jz=interpt（da7（1:41）da7（1:42）t）;
????Jz=Jz*9.8;
????cya=interpt（da1（:1）da1（:2）Ma）;
????mza=interpt（da1（:1）da1（:3）Ma）;
????mzw=interpt（da1（:1）da1（:4）Ma）;
????cx=interpt2（AMCda2aMa）;
else
????Jz=interpt（da7（5:101）da7（5:102）t）;
????Jz=Jz*9.8;
????cya=interpt（da3（:1）da3（:2）Ma）;
????mza=interpt2（TMMda4tMa）;
????mzw=interpt2（TMMda5tMa）;
????cx=interpt2（MCAMMda6aMa）;
end

%一維線性插值
function?b=interpt（xya）;
if（a>=x（1）&&a<=x（end））
???for?i=1:length（x）
???????if（a>=x（i）&&a<=x（i+1））
???????????k=i;
????