用matlab模拟的气枪子波，但是出现问题，求高手
clc;clear;%R2气泡壁质点加速度；
%R1气泡壁质点速度；
%R气泡半径；
n=8;
B=2500;
R(1)=0.084;
R1(1)=0;
y=1.13;
c(1)=1500;
V=0.00246;%(150in3, 1in3=0.00000164m³)
P(1)=13790000;%(2000psi,1psi=6.895kpa)
rou=1025;
dt=1;
R0=0.311;
g=10;
r=1;
h=6;
%p0=P(1).*(R(1)./R0).^(3.*y)
p0=rou*h*g
for t=1:500
P(t+1)=p0.*(R0./R(t)).^(3.*y);
c(t+1)=c(1).*((P(t)+B)./(p0+B)).^((n-1)./(2.*n));
H(t)=(n.*(p0+B).*(((P(t)+B)./(p0+B)).^(1-1./n)-1))./((n-1).*rou);
P1(t)=(-3.*y.*P(t).*R1(t))./R(t);
H1(t)=(((P(t)+B)./(p0+B)).^(-1./n).*P1(t))./rou;
R2(t)=(H(t).*(1+R1(t)./c(t))+R(t).*H1(t).*(1-R1(t)./c(t))./c(t)-1.5.*R1(t).^2.*(1-R1(t)./(3.*c(t))))./(R(t).*(1-R1(t)./c(t)));
dR1(t)=R1(t).*dt+R2(t).*dt.^2./2;
R1(t+1)=R1(t)+dR1(t);
dR(t)=R(t).*dt+R1(t).*dt.^2./2;
R(t+1)=R(t)+dR(t);
w(t)=R(t).*(H(t)+R1(t).^2./2);
K3(t)=c(t).^3.*R(t).^2.*R1(t).*(1-R1(t).^2./(2.*c(t).^2))./w(t).^2+c(t).^2.*R(t).*(1-R1(t)./c(t))./w(t);
u(t)=w(t)./(c(t).*r)+K3(t).*w(t).^2.*(1-w(t)./(c(t).^2.*r)+K3(t).^2.*w(t).^4./(2.*c(t).^8.*r.^4));
p(t)=rou.*(w(t)./r-u(t).^2./2)+(rou.*(w(t)./r-u(t).^2./2).^2)./(2*c(t)^2)+p0;
end
plot(log(p));
