If you call imshow() and directly give it the path to an indexed-color image, it's smart enough to use the colormap from the file instead of gray() or the default axes colormap.
Considering that many of these images I'm seeing on this page are GIF files, that will be what's happening. These are indexed-color images. In order to combine them, you'll convert them to RGB images. Once that's done, there are multiple ways they can be combined.
Let's consider the two images:
You could combine them with imfuse(), but as mentioned, you'll need to convert to RGB first. That said, the only thing it can do is a non-adjustable 50% opacity blend. Its other modes are likely unsuitable for the task of combining color images like these.
[mrpict mrmap] = imread('mr-007.gif');
[tcpict tcmap] = imread('tc-007.gif');
mrpict = ind2rgb(mrpict,mrmap);
tcpict = ind2rgb(tcpict,tcmap);
outpict = imfuse(mrpict,tcpict,'blend');
If alpha compositing suffices, but you just want a different opacity parameter (other than 50%), then you can just do that. Bear in mind that since this is still just a weighted average and the images are largely black, the contributions from the two images will appear darkened.
[mrpict mrmap] = imread('mr-007.gif');
[tcpict tcmap] = imread('tc-007.gif');
mrpict = ind2rgb(mrpict,mrmap);
tcpict = ind2rgb(tcpict,tcmap);
outpict = alpha*tcpict + (1-alpha)*mrpict;
More likely, the task is some form of image blending. Consider the following examples.
[mrpict mrmap] = imread('mr-007.gif');
[tcpict tcmap] = imread('tc-007.gif');
mrpict = ind2rgb(mrpict,mrmap);
tcpict = ind2rgb(tcpict,tcmap);
op2 = 1-(1-tcpict).*(1-mrpict);
op3 = max(tcpict,mrpict);
op4 = sqrt(tcpict) + (1-tcpict).*mrpict;
outpict = [op1 op2; op3 op4];
Take note that the math in the above examples is simplified by the fact that the two RGB arrays are unit-scale floating point. If they were integer-class or if the class could not be known beforehand, those calculations would become more complicated.