function [  ] = leapfrogRAfilter(  )

lambda = 0.8;
h = 1/20;
M = round(4/h);                         % round to the nearest integer
k = lambda*h;                           %
N = round(2.4/k);                       % round to the nearest integer
u0 = setInitialCondition(M);            % Set initial condition

[u1, w1] = LeapFrogRA (M, N, u0, lambda, 0);
[u2, w2] = LeapFrogRA (M, N, u0, lambda, 0.1);
[u3, w3] = LeapFrogRA (M, N, u0, lambda, 0.2);

x = linspace(-1,3,M+1);
plot(x, w1(N-1,:), x, w2(N-1,:), x, w3(N-1,:));
legend('\nu = 0', '\nu = 0.1', '\nu = 0.2');
end


function [u0] = setInitialCondition(M)
% Set an (M+1)-dim row vector containing the initial condition
x = linspace(-1,3,M+1);        % Define (M+1)-dim row vector containing x coordinates
u0 = zeros(1,M+1);             % Allocate memory for row vector u0
ind = (x>=-0.5) & (x <= 0.5);  % Indices of grid points in [-0.5, 0.5]
u0(ind) = cos(pi*x(ind)).^2;   % Set non-zero values for u0
end

function [u, w] = LeapFrogRA(M, N, u0, lambda, nu)
% Implementation of the leap frog scheme with RA filter
u = zeros(N+1, M+1);
u(1,:) = u0;
w = u;
% Use forward-time central space to compute the solution at t_1
u(2, 2:M) = u(1,2:M) - 0.5*lambda*u(1,3:(M+1)) + 0.5*lambda*u(1,1:(M-1));
u(2,1) = 0;
u(2,M+1) = u(2,M);
% The leap frog scheme with RA filter starts from t_2
for n=3:(N+1)
    u(n, 2:M) = w(n-2,2:M) - lambda*( u(n-1,3:(M+1)) - u(n-1,1:(M-1)) );
    u(n,1) = 0;
    u(n,M+1) = u(n,M);
    w(n-1,:) = u(n-1,:) + nu/2*( u(n,:) - 2*u(n-1,:) + w(n-2,:) );
end
end