How to solve the error in creating a panoramic scene from two images? I am using the below code:

S=imread('E:\MADIPALLY\coastal monitoring\project\S_8dec.jpg');
M=imread('E:\MADIPALLY\coastal monitoring\project\M_8dec.jpg');
subplot(1,2,1), imshow(S),title('S')
subplot(1,2,2), imshow(M),title('M')
%%Define the control points (hardcoded for now)
% load controlPoints
fixedPoints = [1589 935;1745 947;1765 510;1288 679];
movingPoints = [989 845;1145 857; 1221 373;677 602 ];
%%Fit geometric transformation to control point pairs
t = fitgeotrans(movingPoints, fixedPoints, 'projective');
%%Apply transformation to moving image and save the spatial location of 
% transformed moving image in the 2-D output space X-Y
[movingTransformed, RBmoving] = imwarp(M, t);
fprintf('The default size for projective transformed image2 is %d X %d X 3 \n',RBmoving.ImageSize)
%%Calculate the XY bounds for the panoramic image

% get the XY bounds of the transformed image in the world coordinates
xData = RBmoving.XWorldLimits;
yData = RBmoving.YWorldLimits;

% get the size of the fixed image
[rowRef,colRef,~] = size(S);

% compute the XY bounds of the panoramic image
xDataOut = [min(1,xData(1)) max(colRef,xData(2))];
yDataOut = [min(1,yData(1)) max(rowRef,yData(2))];
%%Adjust image spatial location for panoramic scene
% Transform both images with the computed panoramic image bounds
% Replace the missing pixels in the image with respect to the panoramic scene, 
% by gray color. R,G and B values for a gray color should be same.

% create a spatial referencing object specifying the panoramic image bounds
% as the XWorldLimits and YWorldLimits. 
Rfixed = imref2d(size(movingTransformed),xDataOut,yDataOut);
% Transform the moving image with the precomputed transform and specifying
% the output image size to be same as the panoramic image size
M = imwarp(M,t,'OutputView',Rfixed,'FillValues',[100 100 100]');

% Transform the fixed image with an identity transform. 
% This would not transform the image but adjust the XY bounds of the image
% to the panoramic image bounds.
t2 = affine2d(eye(3)); % make an identity transform
S = imwarp(S,t2,'OutputView',Rfixed,'FillValues',[100 100 100]');
figure, imshow(S),title('affine transformed image1')
figure, imshow(M),title('projective transformed image2')
fprintf('The computed x/y data size for projective transformed image2 is %d X %d X 3 \n',RBmoving.ImageSize)

%%Create a composite image from the two transformed images
panorama = imfuse(S,M,'blend');
% Display the stitched panoramic scene

figure, imshow(panorama)

M_8dec:

S_8dec: