homography matrix problem
Matrix/vector operations are really strongly optimized in Matlab/Octave. Use them whenever you can.

In your case, instead of multiplying the 3x3 homography matrix with a 3x1 vector for NxM times you can easily modify your calculation to do a multiplication with a 3x3 matrix and a 3xM matrix for N times.

Check out this Octave code:

warning ("off", "Octave:broadcast");
w = 1280;
h = 720;
H = rand(3);

tic();
for i=1:h 
    for j=1:w 
        H*[i;j;1]/(H(3,:)*[i;j;1]); 
    end
end; 
toc();

tic();
for i=1:h 
    v=[ones(1,w)*i;1:w;ones(1,w)];
    (H*v)./(H(3,:)*v); 
end; 
toc();


It gives me the following result:

Elapsed time is 25.5651 seconds.
Elapsed time is 0.273893 seconds.


Maybe you can go further by generating a 3x(N*M) matrix with all the image coordinates.

There is a simple solution to the 'black dots' problem. Right now you have the homography $H$ between the original image and the cover image and you generate the points by looping through the pixels of the original image $(x,y)$ to color the points $(x',y')$ in the destination image.

Let's have the homography $H^{-1}$, loop through all the pixels in the destination image and calculate the pixel address of the original image you want to sample the color from.