Natural boundary condition for 1D heat equation

Hi,

I am trying to solve the 1D heat equation with the following boundary conditions:

I have the follwoing code developed but I need help implementing the natural BC. Can someone please help me?

alpha=2.128*10^-5;
delta_x=0.05;
L=1;
delta_t=10;
t_final=1000;
maxk=t_final/delta_t;
N=L/delta_x;
w=alpha*delta_t/(delta_x^2);
%time grid
t=0:delta_t:t_final;
%mesh grid
x=0:delta_x:L;
[j_min, j_max]=size(x);
% Initial Condition u(x,0)=x/L
for i=1:N+1
    x(i)=(i-1)*delta_x;
    u(i,1)=x(i)/L;
end
% Boundary condition u(0,t)=0, du/dx(L,t)=0
% u(N+1,maxk+1)=zeros;
for Q=1:maxk+1;
    u(1,Q)=0;
end
e = ones(N+1,1);
L1 = spdiags([e -2*e e], -1:1, N+1, N+1);
L1(1,1)=L1(1,1)+1;
L1(N+1,N+1)=L1(N+1,N+1)+1;
% Defining matricies mright and mleft for CN method
aal(1:N-2)=-w;
bbl(1:N-1)=2.+2.*w;
ccl(1:N-2)=-w;
MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);

aar(1:N-2)=w;
bbr(1:N-1)=2.-2.*w;
ccr(1:N-2)=w;
MMr=diag(bbr,0)+diag(aar,-1)+diag(ccr,1);

%implementation of CN method
for P=2:maxk+1    
    uu=u(2:N,P-1);
    u(2:N,P)=inv(MMl)*MMr*uu;
end


figure(1)
plot(x,u(:,end))
xlabel('position')
ylabel('temp')
title('C-N')
%
%
u1(N+1,maxk+1)=zeros;
%Matrix for simple implict 
aa(1:N-2)=-w;
bb(1:N-1)=2+2*w;
cc(1:N-2)=-w;
MM=diag(bb,0)+diag(aa,-1)+diag(cc,1);
%implementation of simple implict 
for B=2:maxk+1;
    for Y=1:maxk;
    u(N+1,Y+1)=u(N+1,Y)+w*(2*u(N,Y)-2*u(N+1,Y));
end
    uuu=u1(2:N,B-1);
    u1(2:N,B)=MM*uuu;
end

figure (2)
plot(x,u1(:,end))
title('Simple Implicit')
xlabel('position')
ylabel('temp')