Calculating Area under Curve

Hello,

I am trying to calculating the area under stress-strain plots by using a code that sums rectangles/triangles for all strain intervals. I know there is a matlab function that does this automatically, such as trapz(), but I need to do it by hand since I am wanting to then decompose into dilational and deviatoric components (Dil + Dev = area under curve). The problem I am running into is I don't know how to code any area above the y-axis as positive and any area under the y-axis as negative. Even though I have strain and stress that gives a linear line, I need this to work for non-liner lines as well. I know this is the problem since I am using the trapz() function to check my code that I developed.

Any direct code to do this would be appreciated, but any general help is also helpful. I have attached my code below.

What I think I need to do is find where my non-linear line crosses my y-axis (where stressxx = 0), however, I don't know how to do this for a non-linear line. I could then write an if loop saying if strainxx greater than where stressx = 0 it would be positive and vice versa.

%clear all; clc; close all
strainxx = [
-1
 2];
stressxx = [
-1
 2];

% Strain Array
for i = 1:length(strainxx)
  % Strain
  strain_array(:,:,i) = [strainxx(i),0,0;0,0,0;0,0,0]; % For Tension
    
  Dil_strain(:,:,i) = 1/3.*trace(strain_array(:,:,i));
  Dil_strain_array(:,:,i) = [Dil_strain(:,:,i), 0, 0; 0, Dil_strain(:,:,i), 0; 0, 0, Dil_strain(:,:,i)];
  Dev_strain_array(:,:,i) = strain_array(:,:,i) - Dil_strain_array(:,:,i);
    
  % Stress
  stress_array(:,:,i) = [stressxx(i),0,0;0,0,0;0,0,0]; % For Tension
  Dil_stress(:,:,i) = 1/3.*trace(stress_array(:,:,i));
  Dil_stress_array(:,:,i) = [Dil_stress(:,:,i), 0,0; 0, Dil_stress(:,:,i), 0; 0,0,Dil_stress(:,:,i)];
  Dev_stress_array(:,:,i) = stress_array(:,:,i) - Dil_stress_array(:,:,i);
end
for n=1:length(strainxx)-1
    
    % Strain
    Dil_delta_strain(:,:,n) = abs(Dil_strain_array(:,:,n+1) - Dil_strain_array(:,:,n));
    Dev_delta_strain(:,:,n) = abs(Dev_strain_array(:,:,n+1) - Dev_strain_array(:,:,n));
    
    Dil_delta_strain1(n) = Dil_delta_strain(1,1,n);
    Dil_delta_strain2(n) = Dil_delta_strain(2,2,n);
    Dil_delta_strain3(n) = Dil_delta_strain(3,3,n);
    
    Dev_delta_strain1(n) = Dev_delta_strain(1,1,n);
    Dev_delta_strain2(n) = Dev_delta_strain(2,2,n);
    Dev_delta_strain3(n) = Dev_delta_strain(3,3,n);
    
    % Stress
    Dil_delta_stress(:,:,n) = abs(Dil_stress_array(:,:,n+1) - Dil_stress_array(:,:,n));
    Dev_delta_stress(:,:,n) = abs(Dev_stress_array(:,:,n+1) - Dev_stress_array(:,:,n));
    
    Dil_delta_stress1(n) = Dil_delta_stress(1,1,n);
    Dil_delta_stress2(n) = Dil_delta_stress(2,2,n);
    Dil_delta_stress3(n) = Dil_delta_stress(3,3,n);
    
    Dev_delta_stress1(n) = Dev_delta_stress(1,1,n);
    Dev_delta_stress2(n) = Dev_delta_stress(2,2,n);
    Dev_delta_stress3(n) = Dev_delta_stress(3,3,n);
    
    
    Dil_stress1(n) = abs(Dil_stress_array(1,1,n));
    Dil_stress2(n) = abs(Dil_stress_array(2,2,n));
    Dil_stress3(n) = abs(Dil_stress_array(3,3,n));
    
    Dev_stress1(n) = abs(Dev_stress_array(1,1,n));
    Dev_stress2(n) = abs(Dev_stress_array(2,2,n));
    Dev_stress3(n) = abs(Dev_stress_array(3,3,n));
    
    Dil_stress1_next(n) = abs(Dil_stress_array(1,1,n+1));
    Dil_stress2_next(n) = abs(Dil_stress_array(2,2,n+1));
    Dil_stress3_next(n) = abs(Dil_stress_array(3,3,n+1));
    
    Dev_stress1_next(n) = abs(Dev_stress_array(1,1,n+1));
    Dev_stress2_next(n) = abs(Dev_stress_array(2,2,n+1));
    Dev_stress3_next(n) = abs(Dev_stress_array(3,3,n+1));
  
    % Slope
    slope(n) = (stressxx(n+1)-stressxx(n))/(strainxx(n+1)-strainxx(n));
    if (stressxx(n) >  stressxx(n+1)) && (stressxx(n) > 0)
        W_dil(n) = 1/2*(Dil_delta_strain1(n)*Dil_delta_stress1(n) + Dil_delta_strain2(n) * Dil_delta_stress2(n) + Dil_delta_strain3(n) * Dil_delta_stress3(n))...
            + Dil_delta_strain1(n)*Dil_stress1_next(n) + Dil_delta_strain2(n) * Dil_stress2_next(n) + Dil_delta_strain3(n) * Dil_stress3_next(n);
       W_dev(n) = 1/2*(Dev_delta_strain1(n)*Dev_delta_stress1(n) + Dev_delta_strain2(n) * Dev_delta_stress2(n) + Dev_delta_strain3(n) * Dev_delta_stress3(n))...
            + Dev_delta_strain1(n)*Dev_stress1_next(n) + Dev_delta_strain2(n) * Dev_stress2_next(n) + Dev_delta_strain3(n) * Dev_stress3_next(n);
     elseif (abs(stressxx(n)) >  abs(stressxx(n+1))) && (stressxx(n) < 0)
        W_dil(n) = 1/2*(Dil_delta_strain1(n)*Dil_delta_stress1(n) + Dil_delta_strain2(n) * Dil_delta_stress2(n) + Dil_delta_strain3(n) * Dil_delta_stress3(n))...
            + Dil_delta_strain1(n)*Dil_stress1_next(n) + Dil_delta_strain2(n) * Dil_stress2_next(n) + Dil_delta_strain3(n) * Dil_stress3_next(n);
       W_dev(n) = 1/2*(Dev_delta_strain1(n)*Dev_delta_stress1(n) + Dev_delta_strain2(n) * Dev_delta_stress2(n) + Dev_delta_strain3(n) * Dev_delta_stress3(n))...
            + Dev_delta_strain1(n)*Dev_stress1_next(n) + Dev_delta_strain2(n) * Dev_stress2_next(n) + Dev_delta_strain3(n) * Dev_stress3_next(n);
    else
        W_dil(n) = 1/2*(Dil_delta_strain1(n)*Dil_delta_stress1(n) + Dil_delta_strain2(n) * Dil_delta_stress2(n) + Dil_delta_strain3(n) * Dil_delta_stress3(n))...
            + Dil_delta_strain1(n)*Dil_stress1(n) + Dil_delta_strain2(n) * Dil_stress2(n) + Dil_delta_strain3(n) * Dil_stress3(n);
       W_dev(n) = 1/2*(Dev_delta_strain1(n)*Dev_delta_stress1(n) + Dev_delta_strain2(n) * Dev_delta_stress2(n) + Dev_delta_strain3(n) * Dev_delta_stress3(n))...
            + Dev_delta_strain1(n)*Dev_stress1(n) + Dev_delta_strain2(n) * Dev_stress2(n) + Dev_delta_strain3(n) * Dev_stress3(n);
    end
    if (slope > 0) && (strainxx(n) < strainxx(n+1)) && (stressxx(n) < stressxx(n+1)) && (stressxx(n) < 0) && (strainxx(n) < 0 ) 
        W_dil(n) = -W_dil(n);
        W_dev(n) = -W_dev(n);
    else (slope > 0) && (strainxx(n) < strainxx(n+1)) && (stressxx(n) < stressxx(n+1)) && (stressxx(n) >= 0) && (strainxx(n) >= 0 )
        W_dil(n) = W_dil(n);
        W_dev(n) = W_dev(n);
    end
end
W_dev_total = sum(W_dev)
W_dil_total = sum(W_dil)
W_total = W_dil_total + W_dev_total
area = trapz(strainxx, stressxx)