nx = 12; ny = 12;
nit=200; 

E = 5;
nu = .3;

A = E/(1 - nu^2);
G = E/(2*(1+nu));


dx = 2/(nx-1); dy = 2/(ny-1);

x=0:dx:2; y=0:dy:2;
[Y,X] = meshgrid(y,x);


u=zeros(nx,ny);un=zeros(nx,ny);
v=zeros(nx,ny);vn=zeros(nx,ny);

C = 1/(2*A/dx^2 + 2*G/dy^2);
D1 = 1/(2*A/dy^2 + 2*G/dx^2);

for iit=1:nit   % iterate for solution
    for i=2:nx-1
        for j=2:ny-1
            % skip the boundary points for inner box
            if i<5 || i>8 || j<5  || j > 8 ||  ...
                (i==5 && j==5)|| (i==5 && j==8)|| (i==8 && j==5)|| (i==8 && j==8)
               un(i,j) = C*(A*(u(i+1,j) + u(i-1,j))/dx^2 ...
                        + G*(u(i,j+1) + u(i,j-1))/dy^2 ...
                        + (G+A*nu)*(v(i+1,j+1) - v(i-1,j+1) ...
                             - v(i+1,j-1) + v(i-1,j-1))/(4*dx*dy));
               vn(i,j) = D1*(A*(v(i,j+1) + v(i,j-1))/dy^2 ...
                        + G*(v(i+1,j) + v(i-1,j))/dx^2 ...
                        + (G+A*nu)*(u(i+1,j+1) - u(i-1,j+1) ...
                        - u(i+1,j-1) + u(i-1,j-1))/(4*dx*dy));
            end
        end
    end

    
    
    % apply downward force across top, y stress = A*dv/dy, so dv/dy = stress/A
    vn(:,ny) = vn(:,ny-1)-1*dy/A;
    
    un(:,ny) = un(:,ny-1);
    
    % natural boundary conditios for sides
    % du/dx =-nu*dv/dy, (u(2,j) - u(1,j))/dx = -nu*(v(1,j+1)-v(1,j))/dy
    
    for k=2:ny-1
        un(1,k) = un(2,k) + nu*(vn(1,k+1)-vn(1,k-1))/2;
        un(nx,k) = un(nx-1,k) - nu*(vn(nx,k+1)-vn(nx,k-1))/2;
        vn(1,k) = v(2,k);
        vn(nx,k) = v(nx-1,k);
        
    end
    
    % boundary conditions for inner box
    
    for k=6:7
        % right edge
        un(8,k) = un(9,k) + nu*(vn(8,k+1)-vn(8,k-1))/2;
        vn(8,k) = vn(9,k);
        
        % left edge
        un(5,k) = un(4,k) - nu*(vn(5,k+1)-vn(5,k-1))/2;        
        vn(5,k) = vn(4,k);
        
        % top edge
        vn(k,8) = vn(k,9) + nu*(un(k+1,8)-un(k-1,8))/2;
        un(k,8) = un(k,9);
        
        %bottom edge    
        vn(k,5) = vn(k,4) - nu*(un(k+1,5)-un(k-1,5))/2;
        un(k,5) = un(k,4);
    end
    
    u=un;
    v=vn;
        
end
                            

u1=u;
v1=v;

um=zeros(nx,ny);
vm=zeros(nx,ny);
for i=1:nx
    for j=1:ny
        um(i,j) = (i-1)*dx+u(i,j);
        vm(i,j) = (j-1)*dy+v(i,j);
    end
end

plot(0,0)
hold on

% plot horizontal lines
for j=1:ny
    if j<=5 || j>=8
       plot(um(:,j),vm(:,j))
    else
        plot(um(1:5,j),vm(1:5,j))
        plot(um(8:12,j),vm(8:12,j))
    end
        
end

%plot vertical lines
for i=1:nx
    if i<=5 || i>=8      
        plot(um(i,:),vm(i,:))
    else
        plot(um(i,1:5),vm(i,1:5))
        plot(um(i,8:12),vm(i,8:12))
    end
end

hold on 

contour(um,vm,shearstress)  % have to class SASplota to calc shearstress
colorbar



        

hold off