drag=0.15; % drag X=[-3.8:0.2:3.8]'; % Initial position Xini=X; Z=interp1(xfull,zfull,X)+rb; U=zeros(size(X)); % Initial velocity EQ=ones(size(X)); EQ(X>xeq2)=2; Xl=size(X,1); col=jet(Xl); nt=890; % number of timesteps dt=0.2; % length of timestep dts=10; % timesteps per save step t=dt*(1:nt); mksz=30; pstm=0.0; pstm1=[]; randon=0; peron=1; randamp=0.2; Arand=randamp*randn(1,nt); peramp=0.4; peromega=0.5; Aper=peramp*sin(peromega*t); Aforc=randon*Arand+peron*Aper; ball_roll_sub