solve pde with variable coefficent

How can i discretize my PDE if it contains a coefficient that include the func. u(x,t) ?

PDE

where e.g.

; IC : and BC

so i tried to use explicit methode until here; assuming

this lee equation from this

https://math.mit.edu/research/highschool/rsi/documents/2017Lee.pdf

but how can i code the a, since i dont know the u(x,t) and solve the pde numerically ?

%% FTCS- ITeration
M=40;             
K=100;            
a=0;             
b=1;              
h=(b-a)/M;       
r=0.49;           
k=r*h^2;          
x=linspace(a,b,M+1);
t=0:k:K*k;       % timestep
U=zeros(K+1,M+1);
uL=0;uR=0; % BC
uIC=sin(4*pi*x); %IC
        U(1,:)=uIC;  U(:,1)=uL; U(:,M+1)=uL; 
        a_u=ones(K+1,M+1); %a is a constant of 1

for kk=2:K+1
    for jj=2:M
        a_p=a_u(kk-1,jj);a_e=a_u(kk-1,jj-1);a_w=a_u(kk-1,jj+1);
        U_W=U(kk-1,jj-1);U_E=U(kk-1,jj+1);U_P=U(kk-1,jj);
        U(kk,jj)= r/4*(a_w-a_e)*(U_W-U_E)+r*a_p*U_W+r*a_p*U_E+(1-2*r*a_p)*U_P;
    end
     U(kk,M+1)=r/4*(a_w-a_e)*(U_W-U_E)+r*a_p*U(kk-1,M)+(1-2*r*a_p)*U(kk-1,M+1);
     U(kk,1)=r/4*(a_w-a_e)*(U_W-U_E)+r*a_p*U(kk-1,2)+(1-2*r*a_p)*U(kk-1,1);
     
end

% plot

   figure
    plot(x,U(1:end,1:end),'linewidth',2);
    a1 = ylabel('Y');
    set(a1,'Fontsize',14);
    a1 = xlabel('x');
    set(a1,'Fontsize',14);
    a1=title(['FTCS Method - r =' num2str(r)]);
    legend('explicite')
    set(a1,'Fontsize',16);
    xlim([0 1]);
    grid;
    
    figure
    [X, T] = meshgrid(x,t);
    s2=surf(X,T,U,'FaceAlpha',0.5,'FaceColor','interp');hold on;
    title(['3-D plot FTCS - r = ' num2str(r)])
    %set(s2,'FaceColor',[1 0 1],'edgecolor',[0 0 0],'LineStyle','--');
    xlabel('x= Length[mm]');ylabel('t= time [s]');zlabel('T = Temperature [°C]');

[Bill Greene]

I suggest using the pdepe function to solve this.

[Muhammad Dzulfikr Islami]

so far i have used this in pdepe func.

it worked, but i just want to code some other method by creating such as BTCS, Crank-nikolson method.

function diff1D

%% h
L = 1; r=0.49;M=40;K=40;
a=0; b=1; h=(b-a)/M; k=4*r*h^2; 
t=0:k:K*k;
x = linspace(0,L,M);

m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
 colormap jet
% imagesc(x,t,sol)
surf(x,t,sol)
% colorbar
xlabel('Distance x','interpreter','latex')
ylabel('Time t','interpreter','latex')
title(['Heat Equation for $0 \le x \le 1$ and $0 \le t $' '  k=' num2str(k)],...
'interpreter','latex')

% u = sol(:,:,1);
function [c,f,s] = heatpde(x,t,u,dudx)
c= 1;
au=(2+u^2);
f= au*dudx;
s=0;
end
    function u0= heatic(x)
        u0=sin(4*pi*x) ;
    end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)

pl = ul;
ql = 0; 
pr = ur;
qr = 0; 
end
end

sd