Numerical solution of an integral equation

Just in the case someone will have the same problem, here is the answer. The problem was in working with vector-argument, vector-value functions in MATLAB, so the kernel should be given in a slightly different way:

function output = kPdfNew(s,t)
z = t-s; 
[nx,ny] = size(s);
output = zeros(size(s));
    for ii = 1:nx
      for jj = 1:ny
        output(ii,jj) = sqrt(2)/pi/(1+(z(ii,jj)+1)^4);
      end
    end
end

Just in case performance and readability are issues, you can achieve what your answer did with the following code as well:

function output = kPdfNew(s,t)
output = sqrt(2)/pi./(1+(t-s+1).^4);
end

(The main problem with your original code was that you iterated for i=1:length(z), instead of length you could have used numel)