Matlab - Crank Nicolson Discretizacion con condiciones de Neumann

hola! amigos, una ayuda, estoy tratando de discretizar este problema con la funcion pico, como un triangulo es la Ecuacion de Difusion 1D con las condiciones de Neumann, no me sale la grafica me sale chueca, e tiene que ir bajando el pico, supuestamente pero me sale como un cuerno! jajaja

tengo dudas en las condiciones de Neumann, alguien me podria decir como las meto? anliticamente ya la hice, pero al pasarlo a la compu ... veo que cambian las cosas

estas son las condiciones de neumann du(0,x)/dx=0 y du(1,x)/dx=0



gracias!

%CRANK NICHOLSON
% Variables, incognitas, pasos y constantes
clc
clear all
%%
n=10 %pasos en x
m=15 %pasos en t
k=0.01 %coeficiente de diffusion
dt=0.1 %paso de tiempo

u=zeros(n+1,m+1) %la matriz donde estaran los datos el recipiente
h=1/n
x=0:h:1
s=k*dt/(2*h^2)

% Datos iniciales y condiciones de borde

%%
%Funcion de Condicion Inicial: Un pico
%se ocupa de todos los valores de la columna 1 Condicion inicial en la matriz u
for i=1:n;
if x(i)<=0.5
u(i)=2.*x(i)
else
u(i)=2-2.*x(i)
end
end

%C.Frontera Neumann
u(n+1,:)=u(n-1,:) %condicion inicial del primer renglon 1 de la matriz u en las J son iguales
u(n,:)=0 %se toman todos los valores de los renglones de la matriz u, en las J son iguales

%matrices del algoritmo
b=ones(n-2,1)
X=eye(n-1)
XX=zeros(1,n-3)
Y=diag(b,-1)

A=-2*eye(n-1)+diag(b,-1)+diag(b,1) %insersion de la diagonal del problema a resolver
% con sus dos diagonales superiores he inferiores
B=inv(eye(n-1)-s*A) %la inversa de A
C=s*[2 zeros(1,n-3) 2]


% Algoritmo del CN
for j=1:m
t(j)=(j-1)*dt %paso del tiempo
u(2:n,j+1)=B*((eye(n-1)+s*A)*u(2:n,j)+C')
end

t(m+1)=m*dt %paso del tiempo
% Grafica
plot(x(1:n+1),u(1:n+1,m))
grid on