Change MATLAB code from Lax-Wendroff to Leapfrog

Question

I want to see how leapfrog would look using this code, but I'm having issues implementing it. I think my biggest problem is adding in the $ U_j^{n-1}$ term, I just don't get the logic.

Here's what the second order leapfrog should look like (correct me if I'm wrong):

$ U_{j+1}^{n} = U_j^{n-1} - \frac{ak}{h}(U_{j+1}^n - U_{j-1}^n) $

Here's the code I'm messing with:

global a
    a = 2;           % advection velocity

    clf              % clear graphics

    ax = 0;
    bx = 1;
    tfinal = 1;                % final time

    h = (bx-ax)/(m+1);         % h = delta x
    k = 0.4*h                  % time step
    nu = a*k/h;                % Courant number
    x = linspace(ax,bx,m+2)';  % note x(1)=0 and x(m+2)=1
                               % With periodic BC's there are m+1 unknowns u(2:m+2)
    I = 2:(m+2);   % indices of unknowns

    nsteps = round(tfinal / k);    % number of time steps
    nplot = 20;       % plot solution every nplot time steps
                      % (set nplot=2 to plot every 2 time steps, etc.)
    nplot = nsteps;  % only plot at final time

    if abs(k*nsteps - tfinal) > 1e-5
       % The last step won't go exactly to tfinal.
       disp(' ')
       disp(sprintf('WARNING *** k does not divide tfinal, k = %9.5e',k))
       disp(' ')
       end

    % initial conditions:
    tn = 0;
    u0 = eta(x);
    u = u0;

    % periodic boundary conditions:
    u(1) = u(m+2);   % copy value from rightmost unknown to ghost cell on left
    u(m+3) = u(2);   % copy value from leftmost unknown to ghost cell on right

    % initial data on fine grid for plotting:
    xfine = linspace(ax,bx,1001);
    ufine = utrue(xfine,0);

    % plot initial data:
    plot(x,u0,'b.-', xfine,ufine,'r')
    axis([0 1 -.2 1.2])
    legend('computed','true')
    title('Initial data at time = 0')

    input('Hit <return> to continue  ');

    % main time-stepping loop:

    for n = 1:nsteps
         tnp = tn + k;   % = t_{n+1}

         % Lax-Wendroff:
         u(I) = u(I) - 0.5*nu*(u(I+1) - u(I-1)) + ...
                       0.5*nu^2 * (u(I-1) - 2*u(I) + u(I+1));

         % periodic boundary conditions:
         u(1) = u(m+2);   % copy value from rightmost unknown to ghost cell on left
         u(m+3) = u(2);   % copy value from leftmost unknown to ghost cell on right

         % plot results at desired times:
         if mod(n,nplot)==0 | n==nsteps
            uint = u(1:m+2);  % points on the interval (drop ghost cell on right)
            ufine = utrue(xfine,tnp);
            plot(x,uint,'b.-', xfine,ufine,'r')
            axis([0 1 -.2 1.2])
            title(sprintf('t = %9.5e  after %4i time steps with %5i grid points',...
                           tnp,n,m+1))
            error = max(abs(uint-utrue(x,tnp)));
            disp(sprintf('at time t = %9.5e  max error =  %9.5e',tnp,error))
            if n<nsteps, input('Hit <return> to continue  '); end;
            end

         tn = tnp;   % for next time step
         end

    %--------------------------------------------------------

    function utrue = utrue(x,t)
    % true solution for comparison
    global a

    % For periodic BC's, we need the periodic extension of eta(x).
    % Map x-a*t back to unit interval:

    xat = rem(x - a*t, 1);
    ineg = find(xat<0);
    xat(ineg) = xat(ineg) + 1;
    utrue = eta(xat);
    return


    %--------------------------------------------------------

    function eta = eta(x)
    % initial data

    beta = 600;
    eta = exp(-beta*(x - 0.5).^2);
    return

Now I believe that it should look something along the lines of....

         % Leap Frog:
         if n==2
             u(I) = utrue(n-1,tnp);
         else
         u(I) = u(I) - nu*(u(I+1) - u(I-1));
         end

But this just isn't working. Also, when run with leapfrog, the code seems to just to not be very accurate(?) as can be seen by the results:

The slopes should align (to show second order accuracy), and the waveform in the first plot should fit the curve more or less.