function?nei
clc;close?all;clear;
i=1;p_=0:0.000003:0.0075;
??S=0.00266;V0=0.00151;m=2.8;??%火炮構造及彈丸諸元
??f=950000;w=1.16;arfa=0.001;roup=1600;theta=0.25;u1=1.7*10^（-8）;n=0.83;
??d=0.00055;
??e1=0.00055;xx=0.75;labda=0.12;uu=0.0;labdas=-0.4104;xxs=1.696;
??p0=3*10^7;%起始條件
??p=p0;
??fai=1.168;%常量計算
??drta=w/V0;
??l0=V0/S;
??ts=0:0.000003:0.0075;
??ffai0=（1/drta-1/roup）/（f/p0+（arfa-1/roup））;
??z0=（（1+4*labda*ffai0/xx）^0.5-1）/2/labda;
??y0=[z0?0?0]‘;
??[ty]=ode45（@funtsy0）;%龍格庫塔法
??Z=y（:1）;
??l=y（:2）;
??v=y（:3）;
???function?dydt=fun（~y）
????????dydt=zeros（31）;??
????????rou=0.2956*（d/2+e1）;
????????zk=（e1+rou）/e1;
???????dydt（2）=y（3）;????????%微分方程2
???????dydt（3）=S*p/fai/m;????%微分方程3
???????if?y（1）????????????dydt（1）=u1/e1*p^n;
???????else
????????????dydt（1）=0;
???????end
????????
???????if?y（1）<1?????
????????????ffai=xx*y（1）*（1+labda*y（