%Set up needed constants
%N=number of spatial intervals
%dx= space step
%dt= timestep
%phi = initial displacement
%psi = initial velocity
%Timesteps = number of Timesteps to take. 
% c= velocitys
%x= vector of spatial discretization points.


N = 20
dx = 0.1
dt = 0.08

c= 1


Timesteps = 30




x=[0:dx:2]

%Define initial data
phi = zeros(1,N+1) %initialize phi, this is redefined below.
psi = zeros(1,N+1)

phi= sin(pi*x)
phi(11:21)=0

s=c^2*dt^2/(dx^2)

%Work matrix contains current values of solution
work = zeros(3,N+1);

work(1,:) = phi

work(2,2:N)=(1-s)*work(1,2:N)+0.5*s*(work(1,1:N-1)+work(1,3:N+1))+2*dt*psi(2:N) 

%diagnostic plot--to be deleteed?
plot(x,work(2,:))
pause


%impose boundary conditions
work(1:2,1) =0 
work(1:2,N+1)=0

for k = 1:(Timesteps-1)
%for j=2:N
%work(3,j)=-work(1,j)+2*work(2,j)+s*(work(2,j+1)+work(2,j-1)-2*work(2,j));
%end
%work(1:2,:) =work(2:3,:);
work(1:2,2:N)=[work(2,2:N);-work(1,2:N)+2*(1-s)*work(2,2:N)+s*(work(2,1:N-1)+work(2,3:N+1))];
plot(x,work(2,:))
axis([0,2,-1.5,1.5])
M(k)=getframe;
end
 
dt*Timesteps

movie(M,1)