Use the patch() function.
Example:
%%% Create sample data
x = 0:1:100;
y = sin(x);
plot(x,y)
patch([x fliplr(x)], [y max(ylim)*ones(size(y))], 'r') %%% Patch the area above the sine curve red
patch([x fliplr(x)], [y min(ylim)*ones(size(y))], 'r') %%% Patch the area below the sine cutve blue
legend() %%% Add the legend here
