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