How can I increase the zero-padding on the DFT?
Your code does not implement the zero padding method. Below is the corrected version that does. Note that zero padding does not improve frequency resolution. It mereley provides interpolated samples of the original spectrum.

   clear all
clc

fs = 50; %50Hz
T = 1/fs;
N = 8;
t = (0:N-1)*T;

y = 5+cos(2*pi*12.5*t)+sin(2*pi*18.75*t);

M = 8;
K = M*N;       % DFT increased size
yk = [y zeros(1,(M-1)*N) ];  %zero padded input

X=zeros(1,K);   % high density

for k = 0:K-1
    for n = 0:K-1
        X(k+1) = X(k+1) + yk(n+1)*exp(-j*(2*pi/K)*k*n);
    end
end

x_mag = abs(X);

f_2 = fs*(0:(K-1))/K;
stem(f_2,x_mag,'-ok','linewidth',2,'MarkerFaceColor','black');