Finite Difference Methods for 1D Transport Equation with Larger Stencil Widths

Question

I'm attempting to re-create the numerical experiments in this paper for solving the 1D scalar advection (transport) equation$:$

$$ \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0 $$

with initial condition $ u(x,0) = e^{-(0.3x)^{2}}$, by using Finite Differences with increasing stencil width. I would like to solve for the equation for $x \in (0, 1050)$ at $t = 1000$.

I have looked through some textbooks on the matter, and the methods I've seen come up most often are explicit methods like the following (taken from the book "Scientific Computing with MATLAB and Octave" by Quarteroni, Saleri, Gervasio). Here, $x_{j} = j \Delta x$, $t^{n} = n \Delta t$, and $\lambda = \Delta t / \Delta x$, and $u_{j}^{n} = u(x_{j}, t^{n})$ are approximate values.

Forward Euler/centered: $u_{j}^{n+1} = u_{j}^{n} - \frac{\lambda}{2}(u_{j+1}^{n} - u_{j-1}^{n})$ which is always unconditionally unstable.

Lax-Friedrichs: $u_{j}^{n+1} = \frac{1}{2}(u_{j+1}^{n} + u_{j-1}^{n}) - \frac{\lambda}{2} (u_{j+1}^{n} - u_{j-1}^{n}) $ which is conditionally stable with respect to a certain norm if the CFL condition is satisfied: $\vert \lambda \vert \leq 1$.

These methods use the Finite Difference approximation with stencil width 3 for the spatial derivative:

$$ f'(x_n) \approx \frac{-\frac{1}{2} f(x_{n-1}) + \frac{1}{2} f(x_{n+1})}{\Delta x} $$ In practice, I wasn't able to get Forward Euler/centered to work, but I was able to get Lax-Friedrichs to work, as in the following implementation in MATLAB. I assumed artificially that we have access to the value of the solution at $x = 0$ for any t, and to deal with the last grid point, I used a linear extrapolation formula:

%%% SOLVING THE 1-D SCALAR TRANSPORT EQUATION %%%
clear all
close all

N = 1049;
x_a = 0;
x_b = 1050;
h = (x_b-x_a)/(N+1);
x = x_a:h:x_b;
dt = 1;
tend = 1000;

lambda = dt/h;

u_t0 = @(x) exp(-(0.3*x).^2);
u_sol = @(x,t) exp(-0.3*(x-t).^2);

u_old = u_t0(x);

time = dt;
while time <= tend
  for i=1:N
    u_new(i+1) = (1/2)*(u_old(i+2) + u_old(i)) - (lambda/2)*(u_old(i+2) - u_old(i));
  endfor
  u_new(1) = u_sol(x_a, time);
  u_new(N+2) = 2*u_new(N-1) - u_new(N);
  time = time + dt;
  u_old = u_new;
endwhile


grid = x_a:0.0001:x_b;
plot(grid, u_sol(grid,tend), '-')
xlim([950 1050])
hold on
plot(x, u_new, '-*')

My question is: how can this be extended to different stencil widths for the Finite Difference formula? For example, if I would like to use the approximation with stencil width 5:

$$ f'(x_n) \approx \frac{\frac{1}{12} f(x_{n-2}) - \frac{2}{3} f(x_{n-1}) + \frac{2}{3} f(x_{n+1}) - \frac{1}{12} f(x_{n+2})}{\Delta x} $$

How could I modify the above methods for this case (or for stencils of arbitrary width)? I imagine this would significantly increase complexity since we would need to approximate the first two and last two grid points by some other method. Or would I have to use a completely different discretization- if so, which one?

Thank you.