Sorry for the late reply ... I just saw this post.
f1 = @(y) g - sign(y).*cd/m.*y.^2;
f2 = @(y) g - sign(y).*(cd/m).*y.^2-(k/m).*(x-L)-(gamma/m).*y;
clearly indicate that the derivatives depend on both the current states of x and y, not just y. So you have a fundamental flaw in the construction of the derivative function handles and things will never work. Since you have two states, x and y, you should pass in both states to your function handles. To be consistent, I would do this for both f1 and f2 even though x is not in f1. E.g.,
f1 = @(x,y) g - sign(y).*cd/m.*y.^2;
f2 = @(x,y) g - sign(y).*(cd/m).*y.^2-(k/m).*(x-L)-(gamma/m).*y;
Having said that, it does lead me to wonder if you have these equations coded correctly. Why does f1 depend only on y but f2 depends on both x and y? You should double check this, and maybe post an image of the differential equations you are solving so we can verify this. I suspect that what you have for this might not be correct.
Your looping would then look something like this, using the current values of x and y (which are x(i) and y(i)) for your derivative calls:
n = numel(t);
x = zeros(size(t));
y = zeros(size(t));
x(1) = x1;
y(1) = y1;
for i = 1:n-1
k1 = h * f1( x(i) , y(i) );
z1 = h * f2( x(i) , y(i) );
k2 = h * f1( x(i) + k1/2, y(i) + z1/2 );
z2 = h * f2( x(i) + k1/2, y(i) + z1/2 );
k3 = h * f1( x(i) + k2/2, y(i) + z2/2 );
z3 = h * f2( x(i) + k2/2, y(i) + z2/2 );
k4 = h * f1( x(i) + k3 , y(i) + z3 );
z4 = h * f2( x(i) + k3 , y(i) + z3 );
x(i+1) = x(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
y(i+1) = y(i) + (z1 + 2*z2 + 2*z3 + z4)/6;
Note that by using spacing I have made the code much more readable and easier to debug.
The variables to plot will be the vectors x and y.
Finally, you should put comments on every line with a constant, stating the units of that constant. This makes it much easier to find unit mismatch errors and helps the reader decipher what the code is doing.