Numerical solution of an EQM for dry friction

Hello,

I have the equation of motion for dry friction case, and I need a solution for it. I impliment Runge_Kutta method to do the solution and write the following .m files but it seems not working!any ideas?

%%Runge-Kutta method for solving equation of motion
% set the initial conditions and the time step
dt=.002;
u0=60/1000;d(1)=u0;
v0=0;v(1)=v0;
% a0=0;
m=2358.9;
k=107000;
zeta= .09;
wn=sqrt(k/m);
c=zeta*(2*m*wn);
wD=wn*sqrt(1-zeta^2);
t=0:dt:5;
 a(1)=accel(v(1),d(1));

for i=1:length(t)-1
    vel=v(i); dis=d(i);
    dv1=dt*accel(vel,dis);
    dx1=dt* vel;
    vel=v(i)+dv1/2; dis=d(i)+dx1/2;
    dv2=dt*accel(vel,dis);
    dx2=dt* vel;
    vel=v(i)+dv2/2; dis=d(i)+dx2/2;
    dv3=dt*accel(vel,dis);
    dx3=dt* vel;
    vel=v(i)+dv3; dis=d(i)+dx3;
    dv4=dt*accel(vel,dis);
    dx4=dt* vel;
    v(i+1)=v(i)+(1/6)*(dv1+2*dv2+2*dv3+dv4);
    d(i+1)=d(i)+(1/6)*(dx1+2*dx2+2*dx3+dx4);
    a(i+1)=accel(v(i+1),d(i+1));
  end

%%For the function accel 
function A = accel(V,D)
m=2358.9;k=107000;zeta= .09;
wn=sqrt(k/m);c=zeta*(2*m*wn);
mue=5;g=386.4;
     A =(-1/m)*(mue*m*g*sign(V)+k*D);
end

Thanks

[Walter Roberson]

What difficulty do you observe? What do you observe that you think should be different?