From the integral that you use in order to calculate marginal PDFs, I guess that x,y are into the interval [0, +inf].
clearvars; clc; close all;
syms x y
Fxy(x,y) = 1-exp(-x)-exp(-y)+exp(-(x+y+x*y)); %Given CDF of X & Y
JPDFxy(x,y) = (x+y+x*y)*exp(-(x+y+x*y)); %Calculated Joint PDF of X & Y
MarginalPDFx(x) = int(JPDFxy,y,0,inf) %Formula for Marginal PDF of X
MarginalPDFy(y) = int(JPDFxy,x,0,inf) %Formula for Marginal PDF of Y
r=0:0.1:10; %plot PDFs for 0<x,y<10
pdfx=zeros(size(r));
pdfy=zeros(size(r));
for n=1:length(r)
pdfx(n)=MarginalPDFx(r(n));
pdfy(n)=MarginalPDFy(r(n));
end
figure;
subplot(2,1,1); plot(r,pdfx,'-b'); grid on; legend('PDFx');
subplot(2,1,2); plot(r,pdfy,'-r'); grid on; legend('PDFy');
Since Fxy is symmetrical, two PDFs are indentical:
