% Usage: [u x t] = wave(f,T,n,nt) or u = wave(f,T,n,nt)
% Wave equation on [0,1] with finite differencing,
% with u(0,t)=u_x(1,t)=0, u(x,0)=f(x), and u_t(x,0)=0
%
% Input:
% f - Matlab inline function f(x) for the initial condition
% T - final time moment
% n - number of subintervals (panels)
% nt - number of time intervals
%
% Output:
% u - computed solution
% x - nodes of the spacial grid
% t - time steps
%
% Examples:
% [u x t] = wave(@smoothbump,2,400,800);          here 'bump' is a user-defined function in M-file
% [u x t] = wave(@smoothbump,2,400,790);
% [u x t] = wave(@smoothbump,2,400,810);

function [u x t] = wave(f,T,n,nt);
x = (0:n)/n;                            % grid on the interval [0,1]
t = (0:nt)*T/nt;                        % grid on the time interval [0,T]
u = zeros(n+1,nt+1);                    % initializing the solution matrix
h = 1/n;                                % spatial mesh-size
tau = T/nt;                             % time step size
c = (tau/h)^2;                          % setting up some constants
c1 = 2-2*c;                             
c2 = c*2;                               
for i=2:n                               
    u(i,1) = f((i-1)*h);                % loop to set up the initial condition at the interior nodes (time step k=1)
end;
u(1,1) = 0;                             % initial condition at the leftmost node x=0 (k=1)
u(1,2) = 0;                             % initial condition at the leftmost node edge x=0 (time step k=2)
u(n+1,1)=f(1);                          % initial condition at the rightmost node x=1 (k=1)
u(n+1,2) = (1-c) * u(n+1,1) + c * u(n,1); % initial condition at the rightmost node x=1 (k=2)
for i=2:n                                                    
    u(i,2) = (1-c) * u(i,1) + c/2 * (u(i+1,1)+u(i-1,1)); % initial condition at the interior nodes (k=2)
end;
for k=2:nt-1
    for i=2:n
        u(i,k+1) = c1 * u(i,k) + c * (u(i+1,k)+u(i-1,k)) - u(i,k-1); % advancing time k -> k+1 at the interior nodes
    end;
    u(1,k+1) = 0;                                                    % advancing time at the leftmost node x=0
    u(n+1,k+1) = c1 * u(n+1,k) + c2 * u(n,k) - u(n+1,k-1);           % advancing time at the rightmost node x=1
end;