function [] = linear_wave()

Nx = 50;
dx = 1/Nx;
x = [0:Nx]'*dx;

c = 1;
T = 1;
Nt = 50;
dt = T/Nt;

u0 = f(x);

% need to solve for u1
u1 = zeros(size(x));
u1(2:end-1) = u0(2:end-1) + ...
    (c*dt/dx)^2/2*(u0(3:end)-2*u0(2:end-1) + u0(1:end-2)) + ...
    dt*g(x(2:end-1));

for nt = 2:Nt
    % centered differencing in time and space
    u = zeros(size(x));

    u(2:end-1) = 2*u1(2:end-1) + ...
    (c*dt/dx)^2*(u1(3:end)-2*u1(2:end-1) + u1(1:end-2)) - u0(2:end-1);
    
    figure(1)
    plot(x,u,'b')
    title(['t = ',num2str(nt*dt)]);
    axis([0 1 -1.5 1.5])
    xlabel('x'); ylabel('u');
    pause(0.1);
    
    u0 = u1;
    u1 = u;
end


%%% Initial Conditions

function [sol] = f(x);
sol = exp(-50*(x-0.5).^2) - exp(-50*0.5^2);


function [sol] = g(x);
sol = zeros(size(x));