Need help coding Finite Difference Scheme for a partial differential equation with both time and space

Hi I need help coding a finite difference scheme for a partial differential equation. My math is correct and I have checked it with my professor but I cannot for the life of me get this code to work . I've included a picture of what the graph is supposed to look like. Also Included my code below.

Question:

Use a finite difference and Runge Kutta scheme to solve the following equation and then code the scheme and graph t=0 and t=2

Equation

Derived Finite difference scheme (Checked with Professor)

Runge Kutta Scheme for Time

Note: n is meant to represent the time space not an exponential

f represents the spatial discretization function

---

Solition Equation

---

Initial Conditions

Code:

clear; clc; close all; 
global Nx Nt dx
%% Initial variables 

tspan = [0 2]; 

dx = 0.1; 
x = -10:dx:10; 

Nx = length(x); 
v = 20; 
x0 = 0; 

CFL = (4*sqrt(2)*sqrt(3))/9; 
fprintf('%d\n',CFL)

dt = dx^3*CFL; 
dt = 0.001; 
t = 0:dt:2; 
Nt = length(t); 

u = zeros(length(x),length(t)); 
u1 = u; 

%% solition solution 

for i = 1:length(t)
    for j = 1:length(x)

u1(j,i) = -v./(2*cosh(1/2.*sqrt(v).*(x(j)-v.*t(i)-x0))^2); 
    end
end

%% initial values case 1 
u0 = u1(:,1); %u1(x,0); 
u(:,1) = u0; 

Loop for Runge Kutta Time space

for j = 1:Nt-1
    
    uj1 = u(:,j); 
    
   
    
    a1 = dt*f5(uj1); 
    a2 = dt*f5(uj1+1/2*a1); 
    a3 = dt*f5(uj1+1/2*a2); 
    a4 = dt*f5(uj1+a3); 
    
    u(:,j+1) = uj1+1/6*(a1+2*(a2+a3)+a4); 
    
     %peridoic boundary 
     u(1,j+1) = u0(1); 
     u(end,j+1) = u0(end); 
    
end

%% Graph
figure()
hold on 
plot(x,u1(:,1))
plot(x,u(:,10),'r')
hold off

Function For Spacial Disc

I used an iterative function but I tried other ones too that I will attach below

function fufu = f5(u)

global Nx dx
u2 = zeros(Nx,1); 

for i = 1:Nx
    
    if (i == 1) %use forward diff
        a1 = (u(i+1)-u(i))/dx; 
        b1 = 1/(2*dx^3)*(-5*u(i)+18*u(i+1)-24*u(i+2)+14*u(i+3)-3*u(i+4)); 
        u2(i) = 6*u(i)*a1 - b1; 
    elseif i == 2 
        a1 = (u(i+1)-u(i-1))/(2*dx); 
        b1 = 1/(2*dx^3)*(-5*u(i)+18*u(i+1)-24*u(i+2)+14*u(i+3)-3*u(i+4)); 
        u2(i) = 6*u(i)*a1 - b1;
    elseif (i == Nx) 
        a1 = (-u(i-1)+u(i))/dx; 
        b1 = 1/(2*dx^3)*(5*u(i)-18*u(i-1)+24*u(i-2)-14*u(i-3)+3*u(i-4)); 
        u2(i) = 6*u(i)*a1 - b1;
    elseif i == Nx-1
        a1 = (u(i+1)-u(i-1))/(2*dx); 
        b1 = 1/(2*dx^3)*(5*u(i)-18*u(i-1)+24*u(i-2)-14*u(i-3)+3*u(i-4)); 
        u2(i) = 6*u(i)*a1 - b1;
    else
        a1 = (u(i+1)-u(i-1))/(2*dx); 
        b1 = (u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2))/(2*dx^3); 
        u2(i) = 6*u(i)*(a1) - b1;
    end
    
     
end

Other Functions I've Tried:

I tried using matrix inversion, other forms of iterative and even a vectorized approach but nothing seems to work

function udot = f(u)

global Nx dx
u0 = u; 

uj = zeros(size(u)); 
for i = 3:Nx-3
    uj(i) = u(i)*6*((u(i+1)-u(i-1))/(2*dx)) - (u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2))/(2*dx^3); 
end
uj(1:2) = u(1:2); 
uj(Nx-1:Nx) = u(Nx-1:Nx); 
uj(1) = uj(end); 

udot  =  uj; 
end

function uFun = f2(u)

global Nx dx

u0 = u; 

a = 1/(2*dx); 
dd(1:Nx) = 0; 
du(1:Nx-1) = 1; 
dl(1:Nx-1) = 1; 
A = diag(dd)+diag(dl,-1)+diag(du,1); 
A = sparse(A); 

b = zeros(Nx,1); 
b(1) = u0(end); 
b(2:end) = u0(2:end)/(6*dx); 

uFun = A\b; 

end


%1st derivative 
%6u*ux = whatever 

function u2 = ufu(u)

global Nx dx

a = 1/(2*dx); 
c = 1/dx^2; 

dd(1:Nx) = 0; 
d1(1:Nx-1) = dx^2*6*u(2:Nx) + 2; 
d2(1:Nx-1) = -dx^2*6*u(2:Nx)- 2;
d3(1:Nx-2) = -1; %(u(i+1))
d4(1:Nx-2) = 1; 

A = diag(dd)+diag(d1,1)+diag(d2,-1)+diag(d3,2)+diag(d4,-2); 
A = sparse(A); 

b = zeros(Nx,1); 
b(1:Nx) = 2*dx^3*u(1:Nx); 



u2 = A\b; 
u2(1) = u2(end); 
end

% Try just the simplified version 
function ufun = f3(u)

global Nx dx

ut = diff(u); %du/dt
ut =[ut;0]; 

b = 1/(2*dx^3); 

dd(1:Nx) = 0; 
dL1(1:Nx-1) = -2*b; 
dL2(1:Nx-2) = b; 
dR1(1:Nx-1) = 2*b; 
dR2(1:Nx-2) = -b; 

A = diag(dd)+diag(dL1,-1)+diag(dL2,-2)+diag(dR1,1)+diag(dR2,2); 
A = sparse(A); 
A(:,1) = 0; 
A(1,1) = 1; 
A(:,end) = 0; 
A(end,end) = 1; 

b1 = zeros(size(u)); 
b1(1:Nx) = ut(1:Nx); 
b1(1) = b(end); 

ufun = A\b1; 

end 

% Trying iterative method again 

function UT = f4(u)
global Nx dx 

dudt = diff(u); 
dudt = [dudt;0]; 
u2 = zeros(Nx,1); 

for i = 1:length(u)
    
    if i == 1
        u2(i) = u(end); 
%         a1 = (u(i+1)-u(end-1))/(2*dx); 
%         u2(i) = 6*u(i)*a1 - (u(i+2)-2*u(i+1)+2*u(end-1)-u(end-2))/(2*dx^3);
    elseif i == 2
        a1 = (u(i+1)-u(i-1))/(2*dx); 
        u2(i) = 6*u(i)*a1 - (u(i+2)-2*u(i+1)+2*u(i-1))/(2*dx^3); 
    elseif i == Nx-1
        a1 = (u(i+1)-u(i-1))/(2*dx); 
        u2(i) = 6*u(i)*a1 - (-2*u(i+1)+2*u(i-1)+u(i-2))/(2*dx^3); 
    elseif i == Nx
        u2(i) = u(end); 
%         a1 = (u(2)-u(i-1))/(2*dx); 
%         u2(i) = 6*u(i)*a1 - (u(3)-2*u(2)+2*u(i-1)+u(i-2))/(2*dx^3);
    elseif (i > 2) && (i < Nx-1)
         a1 = (u(i+1)-u(i-1))/(2*dx); 
         b1 = (u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2))/(2*dx^3); 
         u2(i) = 6*u(i)*(a1) - b1;
    end 
    
    
end 
  
% u2(1) = u2(end); 
  UT = u2;   
    
    

end

What the correct code should look like