Vectorizing Two Sample Kolmogorov-Smirnov test
I'll take this in chronological order, from the top down. I think I have covered all aspects of the code =)

function [OUT] = k_stest(x1,x2,p)


This looks nice, but it would be better if you had spaces between the input arguments. Also, k_stest is not a very good name, as it's hard to understand what it does just be looking at the name. I read in another review that "You don't need the brackets around the output variable, but you should always include it". I don't necessarily agree,  I usually omit them, but that's just a matter of personal preferences.



% Set default p-value
    if nargin == 2 
        p = 0.05;
    end 


This is OK. Depending on how you call the function you might add some additional checks. Is p a scalar? Are x1 and x2 vectors? This is not necessary, but can be useful if you want to use this function a long time from now and don't remember the specifics. assert can be very useful here.



% Compute lenghts
    N1 = length(x1); 
    N2 = length(x2);


numel (number of elements) is better than length. It's more robust and a lot faster for long vectors.



% Set vector t over which the empirical distribution function will be computed.
    t = linspace(min(min([x1 x2])),max(max([x1 x2])),2*max([N1 N2]));


This is dangerous, and a bit messy.

[x1 x2] assumes both x1 and x2 are row-vectors, i.e. horizontal vectors. A lot of functions in Matlab returns vectors as column-vectors i.e. vertical vectors by default, including the colon operator (:). 

If you happen to call this function with two vertical vectors instead of two horizontal vectors, this would fail (due to [x1 x2]).

Instead, you should do:

x12 = [x1(:); x2(:)]; 
t = linspace(min(x12), max(x12), 2*max([N1, N2]));


This way you're converting x1 and x2 to vertical vectors, and concatenate them vertically, prior to the linspace call. Now you only need one call to max and min.



The way you're initializing D is OK!



% Set c from tabulated values
    if p == 0.10
        c = 1.22;
    elseif p == 0.05
        c = 1.36;
    ...


Comparing floating point values using == is dangerous. For instance 0.1 + 0.2 == 0.3 will return false. Instead, you should compare the values with the desired value, but with some tolerance (for instance eps).

    if abs(p - 0.10) &lt; eps
        c = 1.22; 
    elseif abs(p-0.05) &lt; eps
        c = 1.36;


If two values are closer to each other than +/- eps then they are essentially equal.



% Search though the vector t, computing the empirical distribution function for each t, 
% and overwriting the statistic if a higher value is found.
    for i = 1:length(t)
        F1 = sum(x1&lt;=t(i))/N1;
        F2 = sum(x2&lt;=t(i))/N2;
        if abs(F2-F1) &gt; D
            D = abs(F2-F1);
        end
    end


You are right, vectorization is the way to go! Use bsxfun to create a logical matrix where each element that satisfy the condition x1 &lt;= t(ii) are true and the others are false. Then sum the entire mask, and take the maximum difference.

mask1 = bsxfun(@le, x1(:), t);
mask2 = bsxfun(@le, x2(:), t);
s1 = sum(mask1) / N1;
s2 = sum(mask2) / N2;
D = max(abs(s1 - s2));




% Compare the statistic to determine if the samples pass or fail.
    if D == -1;
       disp('Error, invalid input vectors')
       return
    elseif D &gt; c*sqrt((N1+N2)/(N1*N2));
        OUT = 0;
    else
        OUT = 1;
    end    
end


The check for invalid input vectors should be in the start of the function. I can't see a reason why you want it in the end instead. And as far as I can tell it will never fail (if the input is invalid then you will get an error further up).

Instead of the if and elseif, you can simply do:

OUT = D &lt;= c*sqrt((N1+N2)/(N1*N2));




To sum it all up, your function can be written like this. I've omitted the comment for simplicity. You should include them in your code.

function OUT = k_stest(x1, x2, p)

if nargin == 2
    p = 0.05;
end

% Add check for invalid input vectors!

N1 = numel(x1);
N2 = numel(x2);

x12 = [x1(:); x2(:)]; 
t = linspace(min(x12), max(x12), 2*max([N1, N2]));

if abs(p - 0.10) &lt; eps
    c = 1.22; 
elseif abs(p-0.05) &lt; eps
    c = 1.36;
% Continue with the rest of the values
%
end

mask1 = bsxfun(@le, x1(:), t);
mask2 = bsxfun(@le, x2(:), t);
s1 = sum(mask1) / N1;
s2 = sum(mask2) / N2;
D = max(abs(s1 - s2));

OUT = D &lt;= c*sqrt((N1+N2)/(N1*N2));

end