drag=0.15; % drag

X=input('Choose starting position for ball (between -4 and 4; defaults to -3.6): ');
if isempty(X)
  X=[-3.6]'; % Initial position
end
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=250; % number of timesteps
dt=0.2; % length of timestep
dts=10; % timesteps per save step
t=dt*(1:nt);

mksz=30;
pstm=0.05; pstm1=[];

randon=0; peron=1;
randamp=0.2;
Arand=randamp*randn(1,nt);
peramp=0.0; peromega=0.5;
Aper=peramp*sin(peromega*t);
Aforc=randon*Arand+peron*Aper;

ball_roll_sub