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);
    % Create All Nodes Matrix (Interior & Exterior)
    m = zeros(n+2);
    c = vec2mat(w,n); % Vector w rearranged to n x n matrix.
    m(2:n+1,2:n+1) = c; % insert matrix c into middle of m matrix
    % m = symmetric matrix 
    
     h = 10/(n+1);
    
    % Compute dw_dx
    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 % rows
        for j = 1:n+1 % columns
        dw_dx(i,j) = (dx(i,j)+dx(i+1,j))/2;
        end 
    end 
    disp(dw_dx)
    
    % Compute dw_dy
    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 % rows
        for j = 1:n+1 % columns
        dw_dy(i,j) = (dy(i,j)+dy(i,j+1))/2;
        end 
    end 
    disp(dw_dy)
    
    pw_se = zeros(n,n);
    
    for i = 1:16
        se_e1(i) = (1000*10^3*0.3*10^-3*h^2)/(2*(1-0.31^2))*((dw_dx(i))^4/4+(dw_dy(i))^4/4+0.5*(dw_dx(i))^2+0.5*(dw_dy(i))^2);
    end 
    disp(pw_se)