Triple Integration without Symbolic Math Toolbox

Can someone please tell me how to conduct triple integration for the following function 'out' which consists of variables x, y and t. Your help is much appreciated. I know I have to use the triplequad function in MATLAB, but it keeps giving me error.

%defining function components and variables
alpha=4; beta=2; DoF=3; rho=0.7;
cb=1; ab=-0.25; sigmab=0.12; 
cr=1; ar=-0.05; sigmar=0.1;

psi1=tinv(1-exp(-alpha*x),DoF);
psi2=tinv(1-exp(-beta*y),DoF);

comp1=1+(psi1^.2/DoF);
comp2=1+(psi2^.2/DoF);
comp3=(psi1^.2+psi2^.2-2*rho*psi1*psi2)/(DoF*(1-rho^2));

constant=(1/(2*sqrt(1-rho^2)))*alpha*beta*DoF*((gamma(DoF/2)/gamma(DoF/2+1/2))^.2);
copnum=(comp1*comp2)^((DoF+1)/2);
copdenom=(1+(psi1^.2+psi2^.2-2*rho*psi1*psi2)/(DoF*(1-rho^2)))^(1+(DoF/2));

%defining the exponential part for Laplace transform
compA=sqrt((cb*ab)^2+2*(sigmab)^2);
explogA=exp(-(1-t)*compA)*(sigmab^2-(cb*ab+compA))/(sigmab^2-(cb*ab-compA));
Anum=compA+cb*ab+(compA-cb*ab)*explogA;
Adenom=(sigmab^2)*(1-explogA);
At=Anum/Adenom;

compB=sqrt((cr*ar)^2+2*(sigmar)^2);
explogB=exp(-(1-t)*compB)*(sigmar^2-(cr*ar+compB))/(sigmar^2-(cr*ar-compB));
Bnum=compB+cr*ar+(compB-cr*ar)*explogB;
Bdenom=(sigmar^2)*(1-explogB);
Bt=Bnum/Bdenom;

%defining the Laplace transform of exponential function
cop=constant*copnum/copdenom
out=exp(-(At+alpha)*x-(Bt+beta)*y)*cop  %the 'Laplace transform'

Answer 1

[Andrew Newell]

The first thing you need to do is rewrite this as a function with the variables you wish to integrate over as parameters.