function [dep,f,reg] = DiffCentraleElastoPlastique(m,k1,k2,p,DT,Nstep,d1)

dep=zeros(Nstep,1);
vit=zeros(Nstep,1);
f =zeros(Nstep,1);
reg=zeros(Nstep,1);

regime=1;
for i=2:Nstep-1
    switch regime
        case 1; % elastic
            f(i)=k1*dep(i);
        case 2; % augmenting plastic state dep(i)>0
            f(i)=k1*d1+k2*(dep(i)-d1);
        case 3; % elastic unloading
            f(i)=k1*(dep(i)-dR);
        case 4; % augmenting plastic state dep(i)<0
            f(i)=-k1*d1+k2*(dep(i)+d1);
    end
    dep(i+1) = DT^2/m * (p(i)-f(i)) + 2*dep(i)-dep(i-1);
    vit(i+1) = (dep(i+1)-dep(i))/DT;
    reg(i)= regime;
    
    switch regime
        case 1
            if vit(i)>0 && dep(i)>d1;  regime = 2; end
            if vit(i)<0 && dep(i)<-d1; regime = 4; end
        case 2
            if  vit(i)<0; regime = 3; dR=dep(i)-f(i)/k1; end
        case 3
            if vit(i)>0 && dep(i)> d1+k1/(k1-k2)*dR; regime = 2; end
            if vit(i)<0 && dep(i)<-d1+k1/(k1-k2)*dR; regime = 4; end
        case 4
            if  vit(i)>0; regime = 3; dR=dep(i)-f(i)/k1; end
    end
end
