内弹道求解程序

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（