How can i implement the damping coefficient into my codes?

clear all

clc

c = [5 40 200];

m = 20;

k = 20;

h = 0.5;

tn = 15;

y = [1 0];

y01 = y(1);

y02 = y(2);

[t,x,y] = RK2516(@dy1,@dy2,tn,y01,y02,h,c,m,k);

disp([t,x,y]);

----------------------

function dy = dy1(t,y1,y2,c,m,k)

dy = y2;

--------------

function dy = dy2(t,y1,y2,c,m,k)

dy = -(c*y2 - k*y1)/m;

--------------------------

    function [t,y1,y2] = RK2516(dy1,dy2,tn,y01,y02,h,c,m,k)
% dy1 = m-file that evaluates the first ODE
% dy2 = m-file that evaluates the second ODE
% tn = range of initial t value to final t value
% y01 = initial value of the first dependent variable
% y02 = initial value of the second dependent variable
% h = step size

t = [0:h:tn]';

n = length(t);

h2 = h/2;

h3 = h/3;

h6 = h/6;

y1 = y01*ones(n,1);

y2 = y02*ones(n,1);

if t(n) < tn

    t(n+1) = tn;

    n = n + 1;

end

for j = 2:n

    k1 = feval(dy1, t(j-1), y1(j-1), y2(j-1));

    q1 = feval(dy2, t(j-1), y1(j-1), y2(j-1));

    k2 = feval(dy1, t(j-1) + h2, y1(j-1) + h2*k1, y2(j-1) + h2*q1);

    q2 = feval(dy2, t(j-1) + h2, y1(j-1) + h2*k1, y2(j-1) + h2*q1);

    k3 = feval(dy1, t(j-1) + h2, y1(j-1) + h2*k2, y2(j-1) + h2*q2);

    q3 = feval(dy2, t(j-1) + h2, y1(j-1) + h2*k2, y2(j-1) + h2*q2);

    k4 = feval(dy1, t(j-1) + h, y1(j-1) + h*k3, y2(j-1) + h*q3);

    q4 = feval(dy2, t(j-1) + h, y1(j-1) + h*k3, y2(j-1) + h*q3);

    y1(j) = y1(j-1) + h6*(k1 + k4) + h3*(k2 + k3);

    y2(j) = y2(j-1) + h6*(q1 + q4) + h3*(q2 + q3);

end

end