function elastic_wave_simplified( )

rho = 1000;
E = 100;
N = 201;
L = 4;
xp = linspace(0,L,N);

kh = 1;
h = kh*(xp(2)-xp(1));

vp(1:N) = 0;
sp(1:N) = 0;

m = L*rho/(N-1);

t = 0;
dt = 0.01/kh;
tfinal = 20;

while(t<tfinal)
    t = t+dt;

    %calculate spatial derivative of stress
    dsdxp(1:N) = 0;
    for i=1:N
        xi = xp(i);
        si = sp(i);

        dsdx = 0;
        frac = 0;

        for j=1:N
            if (i==j)
                continue;
            end
            xj = xp(j);
            sj = sp(j);
            [~,dwdx] = myzeta(h,xi-xj);
            if (dwdx == 0)
                continue;
            end
            dsdx = dsdx + m/rho*(sj-si)*dwdx;
        end

        dsdxp(i) = dsdx;
    end

    %integrate velocity in time
    vp = vp +dt*1/rho*dsdxp;

    %calculate spatial derivative of velocity
    dvdxp(1:N) = 0;
    for i=1:N
        xi = xp(i);
        vi = vp(i);

        dvdx = 0;
        frac = 0;

        for j=1:N
            if (i==j)
                continue;
            end
            xj = xp(j);
            vj = vp(j);
            [~,dwdx] = myzeta(h,xi-xj);
            if (dwdx == 0)
                continue;
            end
            dvdx = dvdx + m/rho*(vj-vi)*dwdx;
        end

        dvdxp(i) = dvdx;
    end

    sp = sp + dt*E*dvdxp;

    %enforce boundary conditions
    if (t<4)
        sp(N) = 10*sin(pi/4*t);
    else
        sp(N) = 0;
    end

    plot(xp,vp);
    axis([0,4,-0.04,0.04]);
    drawnow;
end

end