i want to solve reaction diffusion,problem is my diffusion coefficient is not constant but a func of spicies itself and it makes the 2 equations coupled i wrote the codes,but im not experienced

as i said it about reaction diffusion equation, i found 2 possible ways for my self: the first one was to use pde toolbox in which i wrote these codes,im not experienced in matlab so please help me out find the problem in my codes, may be there are primary problems even in the definitions so help me out understand them.thanks.

m = 0;
x = linspace(0,1,10);
t = linspace(0,10,2000);
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
%plot(x,u1(end,:))
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [1; 1]; 
f = [0.0933*exp(-5.8*(1-u(2))); 0] .* DuDx; 
F = 2*u(2)*u(1);
s = [0; -F]; 
% --------------------------------------------------------------
function u0 = pdex4ic(x)
u0 = [0; 1]; 
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [0;ul(2)]; 
ql = [1; 0]; 
pr = [0; ur(2)]; 
qr = [1; 0];

the coupled equations are as below: dM/dt=kMH dH/dt=d/dx(DdH/dx) but D itself is: D=D0exp(a(1-M)/1+b(1-M)) with a,b and D0 known constants. D also can be written as D=D0exp(-w(1-M)). the boundary conditions are that at both surfaces there is no diffusion hence surfaces are impermeable.*** i only have this and thats my big problem! i also tried to solve them by difference method as below but i dont know is it even possible or not!

numx = 101;   %number of grid points in x
numt = 2000;  %number of time steps to be iterated
dx = 1/(numx - 1);
dt = 0.00005;

x = 0:dx:1;   %vector of x values, to be used for plotting
t = 0:dt:0.1;   %vector of t values, to be used for plotting

H = zeros(numx,numt);   %initialize everything to zero
D = zeros(numx,numt);
M= zeros(numx,numt);
%specify initial conditions
t(1) = 0;      %t=0     
D0 = 0.000001;
omega= 4.05;
kr=0.194;
W=((0.196-0.025)/2);
U=atanh((0.196+0.025)/(0.025-0.196));
for j=1:numt,i=2:numx-1;
%H(i,1) = exp(-(x(i)-mu)^2/(2*sigma^2)) / sqrt(2*pi*sigma^2);   
H(i,1) = W*tanh((W*3.02*t(j))-U)+W;
D(i,1) = D0*exp(omega*(1-M(i,j)));
M(i,1) = 1-exp(-kr*H(i,j)*t(j));
end

%iterate difference equations
for j=1:numt
t(j+1) = t(j) + dt;
for i=2:numx-1
H(i,j+1) = H(i,j) + (dt/dx^2)*((D(i+1,j)-D(i,j)*(H(i+1,j) - H(i,j))) +D(i,j)*(H(i+1,j)-2*H(i,j)- H(i-1,j))); 
M(i,j+1) = M(i,j)*((kr*H(i,j)*dt)+1);
end
H(1,j+1) = H(2,j+1);          %C(1,j+1) found from no-flux condition
H(numx,j+1) = H(numx-1,j+1);  %C(numx,j+1) found from no-flux condition
end

figure(1);
hold on;
plot(x,H(:,1));
%plot(x,H(:,1001));
%plot(x,H(:,2001));
xlabel('x');
ylabel('H(x,t)');
figure(2);
plot(x,D(:,1));
xlabel('x');
ylabel('D(x,t)');
figure(3);
plot(x,M(:,1));
xlabel('x');
ylabel('M(x,t)');
but answers are really funny. please help me out!