function dh=simulate2d(ip,q,d,deltat);

// transient simulator
// ip and q are row vectors of the same size.
// q is infiltration in meters/day

K=1; W=300; l=300; Sy=0.05; T=50;
// basic hydrolgic data
// T is thickness of aquifer

[M,N]=size(ip);
// input heads as an array M X N
disp(M); disp(N);
dh=zeros(M,N);
for i=1:M // 
  for j=1:N
  flowin=0; 
  if d(i,j)==1 // otherwise do nothing
   
   if (i>1)&(d(i-1,j)>0)
     flowin=flowin+(ip(i-1,j)-ip(i,j))*W*T*K/l;;
     disp(flowin);
   end; 

   if (i<M)&d(i+1,j)>0
     flowin=flowin+(ip(i+1,j)-ip(i,j))*W*T*K/l;
   end;
 
   if (j>1)&d(i,j-1)>0
     flowin=flowin+(ip(i,j-1)-ip(i,j))*W*T*K/l;
   end; 
 
   if (j<N)&d(i,j+1)>0
     flowin=flowin+(ip(i,j+1)-ip(i,j))*W*T*K/l;
   end;
   flowin=flowin+q(i,j)*W*l;
   dh(i,j)=flowin*deltat/(W*l*Sy);
   // the basic transient calculation
  end; // of if
  end; // of inner for
end; // of outer for

return;

