This is quite a tough issue.
Here's a summary of what I've learned so far.
Bio-formats MATLAB plugin's bfsave() can save multiple channels into a *.ome.tif file, which is compatible with ImageJ. But it's quite slow (~2 min per 2048*2048*4 image). Not sure if this is Stitching's problem or bfsave's, but these *.oime.tif files are not compatible with Fiji > plugins > Stitching, as far as I tried.
2. Tiff class
As mentioned in the question, as far as I could see, multi-channel TIFF files saved with Tiff class are not compatible with ImageJ. They are cannot be opened properly. There might be a hidden combination of parameters that allow this to happen. If you find one, let me know!
3. Use imwrite() for each channel
This is a low level work around. Instead of saving four channels at one time, you can save each channel as a separate TIFF and the combine them from within ImageJ.
4. Use ImageJ-MATLAB
As introduced in ImageJ web site, somebody developed a way to commnucate with ImageJ from within MATLAB.
The documentation was incredibly unkind, so I had to spend a lot of time to work out what's going on. Now I've extensively edited the documentation above, so you guys are lucky.
Although it was not very clear in the original documentation, the simple code below can launch an instance of ImageJ from within MATLAB. Moreover, you can throw and catch data between ImageJ and MATLAB bidirectionally, with full access to ImageJ JAVA API.
addpath('/Applications/Fiji.app/scripts') % Update for your ImageJ installation as appropriate
The code below essentially works for me now.
Strangely in my case, the image shown in ImageJ is 90 degrees rotated to the left. Also, the four channels are treated as a stack of four images rather than one image with four channels. So I needed to fix them.
You can also specify LUTs before saving.
This method is ~10 times slower (8.6 sec per 2048*2048*4 image) than saving four separate TIFF files (0.6 sec), but still 10-fold faster than bfsave (~2 min). And you don't have to mess around a folder with many intermediate images.
imp = ij.IJ.getImage();
IJ.run(imp, "Rotate... ", "angle=90 grid=1 interpolation=Bilinear stack");
IJ.run(imp, "16-bit", ""); % convert back to 16 bit per channel
% Need to change stacks to channels
IJ.run("Stack to Images", "");
IJ.run(imp, "Merge Channels...", "c1=-0001 c2=-0002 c3=-0003 c4=-0004 create ignore");
imp2 = IJ.getImage(); % need to get anew
These might be API for stitching.