Some advice, not all directly related to your question, sorry, but I can't help myself
Don't use j as a variable. MATLAB allows j or i as sqrt(-1) but this will eventually bite you when you can't figure out why a script isn't working, and then you realize it's because you redefined sqrt(-1) into something else. (You probably shouldn't use i as a variable either but as we all know j is the true sqrt(-1)!)
Try to define vectors and use ".*" to multiply element-by-element. Much easier then loops.
MATLAB is the only place in the known universe with their weird definition of frequency. Keep sane by using "digital frequency" f = freq_in_Hz / sampling_freq_in_Hz such that 1/2 == Nyquist and follow that convention everywhere, only multiplying by two when you need to use a MATLAB function that thinks 1 == Nyquist. Of course some people prefer omega too. Personally I find it much easier to think in digital [non-angular] frequency and throw in the 2*pi factor wherever I need it.
Now to the stuff that might actually help... To keep things simple, make your LPF with an odd number of taps. Call this L. Then in MATLAB-speak,
M = (L-1)/2; % this is the filter order
n = [-M:M];
hi = h .* cos(2*pi*fc*n); % where fc = desired center frequency / sampling frequency
hq = h .* sin(2*pi*fc*n);
BTW any linear-phase FIR filter of "type 3" or "type 4" automatically has the hilbert-transforming property. That's basically what we're doing when we multiply by sin(2*pi*fc*n), changing from even-symmetric to odd-symmetric at the same time we shift the center frequency away from DC.
Instead of calling fir1 you could do the window design yourself and just specify both filters directly, as
hi = 2*B*sinc(B*n) .* hamming(L)' .* cos(2*pi*fc*n);
hq = 2*B*sinc(B*n) .* hamming(L)' .* sin(2*pi*fc*n);
where L, n, and fc are defined above and B is the desired bandwidth, again in proper normalized units of (true bandwidth / sampling frequency). This is probably the best way to gain insight into what you are doing. (B*sinc(B*n) gives a perfect brick-wall, unity-gain response if n goes on forever, but of course we are truncating it, so we multiply point-by-point with a window function to reduce the sidelobes at the expense of a less-sharp filter, then we modulate from DC to fc by multiplying point-by-point with 2*cos() or 2*sin().
I didn't delve too deeply into your script, but I suspect one of the problems is your LPF bandwidth is too wide and/or you are not modulating the filter far enough in frequency. Remember that multiplication by a real sinusoid produces two copies in the spectral domain, one shifted up and the other shifted down. (Also remember that freqz is only showing you the positive frequencies; the LPF extends the same distance in negative frequency.) If you don't shift far enough, the two copies pile up and interfere with one another. So, the half-band LPF you've designed has a cutoff at 0.25 (I am using proper/correct frequency units here, so 0.25 == halfway to Nyquist). The width of that LPF is 0.5 if you measure its full extent in +/- frequency. When you modulate it you can just barely avoid a spectral collision if you use fc=0.25 (multiply by cos(2*pi*n*0.25)): One copy will be centered at 0.25 and the other at -0.25. This will be essentially an all-pass filter because each copy will span from 0 to 0.5 (or -0.5). Of course, that may be exactly what you want for your SDR.