This Works perfectly @Gianluca, thanks a lot. But the result is still in terms of 'ExpIntegralEi'. Is there a way I can extract the complex function in terms of the variables c1,c2,c3,c4 in a plain form. I tried Series expansion but the function seems to be too complicated.

-((4. (1. E^((1. c1 c3 - 1. c2 c4)/(1. c1 - 1. c2 + 1. c3 - 1. c4))
       ExpIntegralEi[
       0. + 0.5 c2 - 0.5 c3 + (
        0.5 c1 c2 - 0.5 c2^2 - 0.5 c1 c3 + 0.5 c3^2 + 0.5 c2 c4 - 
         0.5 c3 c4)/(c1 - 1. c2 + c3 - 1. c4)] - 
     1. E^((1. c1 c3 - 1. c2 c4)/(c1 - 1. c2 + c3 - 1. c4))
       ExpIntegralEi[
       0. + 0.5 c1 - 0.5 c4 + (
        0.5 c1^2 - 0.5 c1 c2 - 0.5 c1 c3 + 0.5 c2 c4 + 0.5 c3 c4 - 
         0.5 c4^2)/(c1 - 1. c2 + c3 - 1. c4)]))/(
  c1 - 1. c2 + c3 - 1. c4)) + (
 4. (1. E^((1. c1 c3 - 1. c2 c4)/(1. c1 - 1. c2 + 1. c3 - 1. c4))
      ExpIntegralEi[
      0. - 0.5 c2 + 0.5 c3 + (
       0.5 c1 c2 - 0.5 c2^2 - 0.5 c1 c3 + 0.5 c3^2 + 0.5 c2 c4 - 
        0.5 c3 c4)/(c1 - 1. c2 + c3 - 1. c4)] - 
    1. E^((1. c1 c3 - 1. c2 c4)/(c1 - 1. c2 + c3 - 1. c4))
      ExpIntegralEi[
      0. - 0.5 c1 + 0.5 c4 + (
       0.5 c1^2 - 0.5 c1 c2 - 0.5 c1 c3 + 0.5 c2 c4 + 0.5 c3 c4 - 
        0.5 c4^2)/(c1 - 1. c2 + c3 - 1. c4)]))/(
 c1 - 1. c2 + c3 - 1. c4)

Single order approximation

approxintfC = 
     Series[intfC[c1, c2, c3, c4], {c1, 0, 1}, {c2, 0, 1}, {c3, 0, 
       1}, {c4, 0, 1}]

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!!

Thanks for your effort Gianluca, it is a very similar function. I should have checked this myself.

Mathematica ExpintegralEi:

ExpIntegralEi[1.8]
Out[1]= 4.24987

Matlab expint:

>> expint(-1.8)
ans =
  -4.2499 - 3.1416i

but for some reason, - expint(-x) and ExpIntegralEi(x) are matching only in real value and there is an imaginary value associated with the matlab expint that is missing in ExpIntegralEi from Mathematica. But anyway, I shall try out a few other values and see if real(-expint(-x)) does the job. Finally, Thanks a lot Gianluca, you helped me a great deal!!

Does this help?

Assuming[c1>0,
  Integrate[E^(1/4c1(1-y1)(1-y2)),{y1,-1,1},{y2,-1,1}]]