function [seis,p_all] = modeling(nz,nt,source,c,zs,dx,dt,p0,p1)

p = zeros(size(p0));
p_all = zeros(nz,nt);
seis = zeros(nt,1);
seis(1) = p0(zs);
seis(2) = p1(zs);
p_all(:,1) = p0;
p_all(:,2) = p1;
for it=3:nt
    for iz=2:nz-1
        alpha = c(iz)*dt/dx;
        alpha = alpha*alpha;
        p(iz) = 2*(1-alpha)*p1(iz) + alpha*(p1(iz+1)+p1(iz-1)) - p0(iz);
    end;
    p(zs) = p(zs) + source(it);
    seis(it) = p(zs);
    p0 = p1;
    p1 = p;
    p_all(:,it) = p;
end;
