  clear all; clc
  nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl;  
  ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11;
  nt=200; nout=10; h=0.0005; t=0;
  w=ones(nx,ny);
  for i1=1:nout
    for i2=1:nt
      for i=1:nx
        for j=1:ny  
          if(i== 1)      w( 1,j)=0; wt( 1,j)=0;   
            elseif(i==nx)w(nx,j)=0; wt(nx,j)=0; 
            elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
                                 +2*(w(i,2)-w(i,1))/dy^2;
            elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
                         wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
                                 +(wf-2*w(i,ny)+w(i,ny-1))/dy^2;
            else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
                        +(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
          end             
        end 
      end 
      w=w+wt*h;   
      t=t+h;
    end
      fprintf('t = %5.2f   w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
  end
      fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));
      
