Duplicate Equation in MatLab ODE45 Solve

I'm trying to write a MatLab code which can extend the number of ODE's to solve for simply. This is the code that I currently have, for simplicity I am starting with a spring mass damper system.

    clear;clc;close all;
    tspan=[0:0.01:1];
    x0=[1;0;-1;0;7;0;-7;0;5;0;-5;0;10;0;-10;0];
    [t,x] = ode45(@Spring_Mass_Damper,tspan,x0);
    figure(1)
    plot(t,x(:,1));
    hold on
    plot(t,x(:,3));
    hold on
    plot(t,x(:,5));
    hold on
    plot(t,x(:,7));
    hold on
    plot(t,x(:,9));
    hold on
    plot(t,x(:,11));
    hold on
    plot(t,x(:,13));
    hold on
    plot(t,x(:,15));
    grid on
    xlabel('Time')
    ylabel('Displacement(x)')
    title('Displacement Vs Time')

function xp = Spring_Mass_Damper(t,x)
c1=10;
G=9.81;
m1=1;
k1=2000;


xp=[[x(2)];(-(c1/m1).*[x(2)]-(k1/m1).*[x(1)])-(G.*m1);
    [x(4)];(-(c1/m1).*[x(4)]-(k1/m1).*[x(3)])-(G.*m1);
    [x(6)];(-(c1/m1).*[x(6)]-(k1/m1).*[x(5)])-(G.*m1);
    [x(8)];(-(c1/m1).*[x(8)]-(k1/m1).*[x(7)])-(G.*m1);
    [x(10)];(-(c1/m1).*[x(10)]-(k1/m1).*[x(9)])-(G.*m1);
    [x(12)];(-(c1/m1).*[x(12)]-(k1/m1).*[x(11)])-(G.*m1);
    [x(14)];(-(c1/m1).*[x(14)]-(k1/m1).*[x(13)])-(G.*m1);
    [x(16)];(-(c1/m1).*[x(16)]-(k1/m1).*[x(15)])-(G.*m1)];
end

The Main question I have regards this area of code:

xp=[[x(2)];(-(c1/m1).*[x(2)]-(k1/m1).*[x(1)])-(G.*m1);
    [x(4)];(-(c1/m1).*[x(4)]-(k1/m1).*[x(3)])-(G.*m1);
    [x(6)];(-(c1/m1).*[x(6)]-(k1/m1).*[x(5)])-(G.*m1);
    [x(8)];(-(c1/m1).*[x(8)]-(k1/m1).*[x(7)])-(G.*m1);
    [x(10)];(-(c1/m1).*[x(10)]-(k1/m1).*[x(9)])-(G.*m1);
    [x(12)];(-(c1/m1).*[x(12)]-(k1/m1).*[x(11)])-(G.*m1);
    [x(14)];(-(c1/m1).*[x(14)]-(k1/m1).*[x(13)])-(G.*m1);
    [x(16)];(-(c1/m1).*[x(16)]-(k1/m1).*[x(15)])-(G.*m1)];

is there a simple way in which i can achieve the same number of systems(ODE's) without having to copy and paste the first equation in the block and manually change the x indexing?

This should also work:

xp = zeros(size(x,1),1)
for idx=1:size(x,1)
    if mod(idx,2) == 0
        xp(idx) = (-(c1/m1).*[x(idx)]-(k1/m1).*[x(idx-1)])-(G.*m1);
    else
        xp(idx) = x(idx);
    end
end

The plotting routine can be simplified in a similar manner, like this:

figure(1)
hold on
for idx=1:2:size(x,2)
    plot(t,x(:,idx));
end
grid on
xlabel('Time')
ylabel('Displacement(x)')
title('Displacement Vs Time')