%  This code solves the wave equation using centered differences
%  in both space and time.

xmin = 0; xmax = 1;
n = input(' Enter number of spatial subintervals: ');
h = (xmax-xmin)/n;

%  Set initial value.
ukm1 = zeros(n-1,1);
for i=1:n-1,
  xi = i*h;
  ukm1(i) = xi*(1-xi);
end;
usave(:,1) = ukm1;    % Save u at each timestep in array usave.

tmin = 0; tmax = 2;
dt = input(' Enter delta t: ');

asq = 2;  % Set parameter a^2.
c = asq*dt^2/h^2;

%  Use forward Euler for the first time step.
uk = ukm1;
plot([0:h:1]', [0; uk; 0]), axis('equal'), axis([0 1 -.25 .25]), pause(1),
usave(:,2) = uk;

%  Loop over time steps.
k = 1;
ukp1 = zeros(n-1,1);
for t=2*dt:dt:tmax,
  k = k+1;
  ukp1(1) = 2*uk(1) - ukm1(1) + c*(uk(2) - 2*uk(1));
  for i=2:n-2,
    ukp1(i) = 2*uk(i) - ukm1(i) + c*(uk(i+1) - 2*uk(i) + uk(i-1));
  end;
  ukp1(n-1) = 2*uk(n-1) - ukm1(n-1) + c*(-2*uk(n-1) + uk(n-2));
  plot([0:h:1]', [0; ukp1; 0]), axis('equal'), axis([0 1 -.25 .25]), pause(1),
  usave(:,k) = ukp1;
  ukm1 = uk;
  uk = ukp1;
end;