function [pot_en,pot_en1] = energy3d(x,y,z,kappa,lambda)
[m,n] = size(x);
pot_en = 0; pot_en1 = 0;
for i=1:m-1
    for j=1:n-1
            A = [x(i,j),y(i,j),kappa/2*z(i,j)^2-lambda*z(i,j)];
            B = [x(i+1,j),y(i+1,j),kappa/2*z(i+1,j)^2-lambda*z(i+1,j)];
            C = [x(i,j+1),y(i,j+1),kappa/2*z(i,j+1)^2-lambda*z(i,j+1)];
            Z = [A(3),B(3),C(3)];
            AB = B - A;
            AC = C - A;
            N = cross(AB,AC);
            A1 = [x(i,j),y(i,j),z(i,j)];
            B1 = [x(i+1,j),y(i+1,j),z(i+1,j)];
            C1 = [x(i,j+1),y(i,j+1),z(i,j+1)];
            AB1 = B1 - A1;
            AC1 = C1 - A1;
            N1 = cross(AB1,AC1);
            surf_area = sqrt(N1(1)^2+N1(2)^2+N1(3)^2);
            zmin = min(Z);
            zmax = max(Z);
            pot_en = pot_en + surf_area/2 + (sum(Z))/3*abs(N1(3))/2;
            pot_en1 = pot_en1 + surf_area/2 + (2*zmin+zmax)/3*abs(N1(3))/2;
    end
end
for i=2:m
    for j=2:n
            A = [x(i,j),y(i,j),kappa/2*z(i,j)^2-lambda*z(i,j)];
            B = [x(i-1,j),y(i-1,j),kappa/2*z(i-1,j)^2-lambda*z(i-1,j)];
            C = [x(i,j-1),y(i,j-1),kappa/2*z(i,j-1)^2-lambda*z(i,j-1)];
            Z = [A(3),B(3),C(3)];
            AB = B - A;
            AC = C - A;
            N = cross(AB,AC);
            A1 = [x(i,j),y(i,j),z(i,j)];
            B1 = [x(i-1,j),y(i-1,j),z(i-1,j)];
            C1 = [x(i,j-1),y(i,j-1),z(i,j-1)];
            AB1 = B1 - A1;
            AC1 = C1 - A1;
            N1 = cross(AB1,AC1);
            surf_area = sqrt(N1(1)^2+N1(2)^2+N1(3)^2);
            zmin = min(Z);
            zmax = max(Z);
            pot_en = pot_en + surf_area/2 + (sum(Z))/3*abs(N1(3))/2;
            pot_en1 = pot_en1 + surf_area/2 + (2*zmin+zmax)/3*abs(N1(3))/2;
    end
end