内弹道matlab程序

clear;
close all;

format long
d=0.0127;
S=0.82*0.0127^2;
V0=2.04e-5;
l_0=V0/S;
lg=0.924;
f=1000000;
alpha=0.001;
w=0.017;
rou=1600;
theta=0.2;
phi=1.45;
chi=0.79825;
lamda=0.1387;
mu=-0.043956;
e1=0.00052/2;
u1=7.5991e-10;
Is=e1/u1;
chi_s=1.2645;
lamda_s=-0.31322;
zk=1.4434;
%Ik=447000;
m=0.048;
p0=30e6;
delta=800;
psi0=（1/delta-1/rou）/（f/p0+alpha-1/rou）;
sigma0=sqrt（1+4*lamda*psi0/chi）; 
z0=2*psi0/chi/（sigma0+1）;%（sigma0-1）/2/lamda;

%====赋予初值====%?
v（1）=0;
l（1）=0;
p（1）=p0;
z（1）=z0;
psi（1）=psi0;
lpsi（1）=l_0*（1-delta/rou-（alpha-1/rou）*delta*psi（1））;
t（1）=0;
h=0.000001;

for i=1:100000
z1=p（i）/Is;
v1=S*p（i）/m/phi;
l1=v（i）;
psi1=（z（i）>=0&z（i）<=1）.*（chi+2*chi*z（i）*lamda+3*chi*mu*z（i）^2）*z1+（z（i）>1&z（i）zk）.*0;
lpsi1=-l_0*（alpha-1/rou）*delta*psi1; 
p1=（（f*w/S+p（i）*l_0*delta*（alpha-1/rou））*psi1-theta*phi*m*v1*v（i）/S-p（i）*l1）/（l（i）+lpsi（i））;

z2=（p（i）+h*p1/2）/Is;
v2=S*（p（i）+h*p1/2）/m/phi;
l2=v（i）+h*v1/2; 
psi2=（z（i）>=0&z（i）<=1）.*（chi*z2+2*chi*（z（i）+h*z1/2）*lamda*z2+3*chi*mu*z2*（z（i）+h*z1/2）^2）+（z（i）>1&z（i）<=zk）.*（chi_s*z2+2*chi_s*z2*lamda_s*（z（i）+h*z1/2））+（z（i）>zk）.*0;
lpsi2=-l_0*（alpha-1/rou）*delta*psi2;
p2=（（f*w/S+（p（i）+h*p1/2）*l_0*delta*（alpha-1/rou））*psi2-theta*phi*m*v2*（v（i）+h*v1/2）/S-（p（i）+h*p1/2）*l2）/（（l（i）+h*l1/2）+（lpsi（i）+h*lpsi1/2））;

z3=（p（i）+h*p2/2）/Is;
v3=S*（p（i）+h*p2/2）/m/phi;
l3=v（i）+h*v2/2; 
psi3=（z（i）>=0&z（i）<=1）.*（chi*z3+2*chi*（z（i）+h*z2/2）*lamda*z3+3*chi*mu*z3*（z（i）+h*z2/2）^2）+（z（i）>1&z（i）zk