% solve u_xx = c^2 * u_tt
% x := space coordinate [0,a].
% t := time [0,T].
% u := the function to be solved.

% initial condition   : u(x,t=0) = sin(pi*x).
% initial velocity    : u_t(x,t=0) = V(x) = 0.
% boundary conditions : u(x=0,t) = u(x=a,t) = 0.


c = 3.0;
a = 5.0;
T = 10.0;

time_step = 1.0e-3; 
x_step = 1.0e-2;

const = c*time_step/x_step;
const = const*const;
const_half = 0.5*const;
const_1 = 1.0-const;
const2 = 2.0 * const_1;

Ntime = floor(T/time_step) + 1;
Nx = floor(a/x_step) + 1;

%%x = zeros(Ntime,Nx);
x = zeros(Nx,Ntime);

tic;

% set up initial condition:
for jt = 1 : Ntime
%    u(jt,1)  = 0.0;
%    u(jt,Nx) = 0.0;
    u(1,jt)  = 0.0;
    u(Nx,jt) = 0.0;
end

for ix = 2 : Nx-1
%    u(1,ix) = sin(pi*x_step);
%    u(2,ix) = const_half * ( u(1,ix+1) + u(1,ix-1) ) + const_1 * u(1,ix); 
    u(ix,1) = sin(pi*x_step*ix);
    u(ix,2) = const_half * ( u(ix+1,1) + u(ix-1,1) ) + const_1 * u(ix,1); 
end

for jt = 2 : Ntime-1
for ix = 2 : Nx-1
%    u(jt+1,ix) = const2*u(jt,ix) + const*(u(jt,ix+1) + u(jt,ix-1)) - u(jt-1,ix);
    u(ix,jt+1) = const2*u(ix,jt) + const*(u(ix+1,jt) + u(ix-1,jt)) - u(ix,jt-1);
end
end

toc;
