%              wavemodel.m
%              need courant number cn, timesteps nt
u=zeros(1,51);
h=ones(1,50);
a=2:length(h)+1;
b=2:length(h);
freq=4*2*pi/nt
hold off
axis([0,length(a),-.5,1.5])
plot(b,u(b),b,h(b))
hold on
for t=1:nt
     h(1)=1+0.1/cosh(freq*t-1)^2;
     u(b)=u(b)-(u(b).*(u(b+1)-u(b-1))/2+diff(h))*cn;
     h(b)=h(b)-diff(u(a)).*(h(b-1)+h(b))/2;
     if floor(t/10)*10==t
      plot(b,u(b),b,h(b))
     end
end
