nx = 10; ny = 10;
nit=500; 

E = 3;
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
            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
    
    
    %show(un)
    %show(vn)
    %keyboard

    
    % natural boundary conditios for top
    % tauyx = 1
    % G(du/dy+dv/x) = 1, du/dy = 1/G-dv/dx
    % sigmayy = 0
    % dv/dy = -nu*du/dx
    
    
    %vn(1,10) = vn(1,9) - nu*(un(2,10)-un(1,10));
    %un(1,10) = un(1,9) + 1/G - (un(2,10)-un(1,10));
    %vn(10,10) = vn(10,9);
    %un(10,10) = un(10,9);
    %for k=2:9
        %un(k,10) = un(k,9) + 1/G - (un(k+1,10)-un(k-1,10))/2;
        %vn(k,10) = vn(k,9) - nu*(un(k+1,10)-un(k-1,10))/2;     
    %end
    
    % take 2
   
  
    
   
    

    
    % natural boundary conditios for sides
    
    % sigmxx = 0
    % du/dx =-nu*dv/dy, (u(2,j) - u(1,j))/dx = -nu*(v(1,j+1)-v(1,j))/dy
    % tauxy = 0
    % dv/dx = -du/dy, (v(10,j) - v(9,j))/dx = -(u(10,j+1) - u(10,j-1))/2dx
    %for k=2:9
    %   un(1,k) = un(2,k) + nu*(vn(1,k+1)-vn(1,k-1))/2;
    %   un(10,k) = un(9,k) - nu*(vn(10,k+1)-vn(10,k-1))/2;
       
    %end
    un(1,2:9) = un(2,2:9);
    un(10,2:9) = un(9,2:9);
    vn(1,2:9) = vn(2,2:9);
    vn(10,2:9) = vn(9,2:9);
    
    % boundary conditions along top
    vn(:,10) = vn(:,9);
    un(:,10) = un(:,9) + .1;
     
    
    u=un;
    v=vn;
        
end
                            


um=zeros(nx,ny);   % array x coordinates of stressed grid points
vm=zeros(nx,ny);   % array y coordinates of stressed grid points

% calculate the coordinates of the stressed grid points
for i=1:nx
    for j=1:ny
        um(i,j) = (i-1)*dx+u(i,j); % original x coord. + x displacement
        vm(i,j) = (j-1)*dy+v(i,j); % original y coord. + y displacement
    end
end

figure(1)
plot(0,0)
hold on

%plot the horizontal lines between grid points
for j=1:ny
    plot(um(:,j),vm(:,j))  % plot jth horizontal line
end

%plot the vertical lines between grid points
for i=1:nx
    plot(um(i,:),vm(i,:))  % plot ith vertical line
end

hold off