SE_diff.asv

function [pw_se,w] = SE_diff(T,P,n)
    
    % Call function from part c to obtain vector w
    w = membrane_solution(T,P,n);
    m = zeros(n+2);
    c = vec2mat(w,n);
    m(2:n+1,2:n+1) = c;
    
    
    h = 10/(n+1);
    
    w_bar = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
            w_bar(i,j) = (m(i,j)+m(i,j+1)+m(i+1,j)+m(i+1,j+1))/4;
        end
    end
    
    left_side = zeros(n+1,n+1);
    for  i = 1:(n+1)^2
        left_side(i) = P * h^2*w_bar(i);
    end
    L = sum(sum(left_side));
     
    dx = zeros(n+2,n+1);
    for i = 2:n+2
        dx(:,i-1)= (m(:,i)-m(:,i-1))/(h);
    end
 
    dw_dx = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
        dw_dx(i,j) = (dx(i,j)+dx(i+1,j))/2;
        end 
    end 

    dy = zeros(n+1,n+2);
    for i = 2:n+2
        dy(i-1,:)= (m(i-1,:)-m(i,:))/(h);
    end
    
    dw_dy = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
        dw_dy(i,j) = (dy(i,j)+dy(i,j+1))/2;
        end 
    end 

    right_side = zeros(n+1,n+1);
    
    for i = 1:(n+1)^2
        right_side(i) = (0.25*(dw_dx(i))^4+0.25*(dw_dy(i))^4+0.5*(dw_dx(i)*dw_dy(i))^2);
    end 
    R = sum(sum(right_side));
    R = (1000*10^3*0.3*10^-3*h^2)/(2*(1-0.31^2))*R;
    
    pw_se = R - L;
end