L=100;
i=zeros(1,L);
v=zeros(1,L;
dz=0.01;
c=0.1;
l=250;
dt=dz*sqrt(l*c);
for n=1:2000
t=n*dt;
for k=1:L-1
if k==1
v(k)=v(k)-dt/(c*dz)*(i(k+1)-i(k));
elseif k==L-1
v(k)=v(k)-dt/(c*dz)*(i(k+1)-i(k));
else
v(k)=v(k)-(((6/7)*dt)/(c*24*dz))*(i(k-1)-(27*i(k))+(27*i(k+1))-i(k+2));
end
end
v(1)=sin(2*pi*t);
v(st)=0;
for k=2:L
if k==2
i(k)=i(k)-dt/(l*dz)*(v(k)-v(k-1));
elseif k==st
i(k)=i(k)-dt/(l*dz)*(v(k)-v(k-1));
else
i(k)=i(k)-(((6/7)*dt)/(l*24*dz))*(v(k-2)-(27*v(k-1))+(27*v(k))-v(k+1));
end
end
plot(v);
axis([1 L -3 3]);
frame = getframe;
end