I checked right now the Help Center of Matlab. Among the special functions it offers also expint, described as "Exponential integral function". I wonder if it is the same thing as Mathematica's ExpIntegralEi.

I would try this way:

c = {c1, c2, c3, c4};
sf = 1/4*{(1 - y1)*(1 - y2),
    (1 + y1)*(1 - y2),
    (1 + y1)*(1 + y2),
    (1 - y1)*(1 + y2)};
cGP = sf . c;
fC = Exp[cGP];
firstIteratedIntegral = Integrate[fC, {y1, -1, 1}];
primitiveOfFirstIteratedIntegral[y2_] = 
  Integrate[firstIteratedIntegral, y2];
doubleIntegral[c1_, c2_, c3_, c4_] = 
 primitiveOfFirstIteratedIntegral[1] - 
  primitiveOfFirstIteratedIntegral[-1]
Manipulate[Quiet@Plot3D[doubleIntegral[c1, c2, c3, c4],
   {c1, -2, 2}, {c2, -2, 2}, AxesLabel -> Automatic],
 {{c3, -1}, -2, 2},
 {{c4, 1/2}, -2, 2}]

I use Matlab to develop a Finite Element model of a mechanical process. This process is governed by a series of PDEs. The whole volume is broken into small elements and these PDEs are kinda converted into a form called weak form which is integrated over the volume of these elements to calculate a residual which is equated to '0', this is basiclly what FEM is. Usually the integral over the element is calculated through Gaussian Quadrature, But I'm trying to do it using Matematica analytical integration as I'm using the concept of reduced integration for certain terms of the PDE as it demands this. So anyway, this Integrated expression goes into the residual formulation in Matlab as a function of the field variables 'c'. SoI'll need an expression that I can write in matlab to compute the value of the residual at each Newton Iteration for every element. Thanks a lot for your help and curiosity, I hope I gave a clear answer.

TL;DR: I need to write this expression in matlab where I give the c values and get the expression value as an output. Thanks!!