function [hh,x,nn,deltaj]=wellsim(q,L,m,K,Sy,dt,n,h)

// H is the lake height
// q is the rainfall, L is the length of the land
// m is the number of cells
// K is conductance, Sy is the specific yield
// n is the number of time steps, dt is the time step
// eps is the error bound, h is the starting heads.
// hh: output, x:the X axis points, nn:# iter, deltaj: errors

hh=h; ii=1;
delta=1;

dd=L/m;
for i=1:m
  x(i)=L/2/m+(i-1)*dd;
end;

while ii<n

for j=2:m-1
  vin=(h(j+1)-h(j))*K*h(j+1)/dd;
  vin=vin+(h(j-1)-h(j))*K*h(j)/dd;
  deltaj(j)=vin/Sy/dd;
  hn(j)=h(j)+vin*dt/Sy/dd;
end;
vin=(h(m-1)-h(m))*K*h(m)/dd;
hn(m)=h(m)-q*dt/Sy/dd+vin*dt/Sy/dd;
vin=(h(2)-h(1))*K*h(2)/dd;
hn(1)=h(1)+vin*dt/Sy/dd;
h=hn;
ii=ii+1;

end;

hh=h;
nn=ii;

endfunction;

 

