Coeficientes de difusão variável no matlab

Question

Estou ficando sem tempo para esse código, então qualquer ajuda seria muito apreciada. Atualmente, estou codificando a equação de calor / difusão 2D em matlab. minha equação é assim:

∂u∂t=αx∂2u∂x2+αy∂2u∂y2+Au(1−uK)
∂u∂t=αx∂2u∂x2+αy∂2u∂y2+Au(1−uK)

Estou usando o Neumann bc's com o fluxo = 0 em todos os lugares. u é o número da população, A, a taxa de crescimento intrínseca e K a capacidade de carga. O código que estou usando é um código que encontrei aqui e que editei para usar o método implícito / BTCS.

Eu geramos resultados usando valores constantes para os coeficientes de difusão xey. agora preciso alterar o coeficiente de difusão em x para uma função em x e o coeficiente de difusão em y para uma função em y. alguma idéia de como eu posso fazer isso e quais outras seções no código precisam ser corrigidas? Eu já comecei a fazer alterações, deixe-me saber o que mais precisa ser feito.

    %%
%Specifying parameters
nx=50;                           %Number of steps in space(x)
ny=50;                           %Number of steps in space(y)       
nt=500;                          %Number of time steps 
dt=0.01;                         %Width of each time step
dx=2/(nx-1);                     %Width of space step(x)
dy=2/(ny-1);                     %Width of space step(y)
x=0:dx:2;                        %Range of x(0,2) and specifying the grid points
y=0:dy:2;                        %Range of y(0,2) and specifying the grid points
u=zeros(nx,ny);                  %Preallocating u
un=zeros(nx,ny);                 %Preallocating un
visxstr = input('x diffusion coefficient ','s');
visystr = input('y diffusion coefficient ','s');
visx = inline(visxstr);          %Diffusion coefficient/viscocity
visy = inline(visystr);          %Diffusion coefficient/viscocity
RepRateA = 0.3;                  %Intrinsic Growth Rate of Population
CarCapK = 1.7;                   %Carrying Capacity of Environment 
f = RepRateA*u(1-u/CarCapK);     %Nonlinear Source Term
UnW=0;                           %x=0 Neumann B.C (du/dn=UnW)
UnE=0;                           %x=L Neumann B.C (du/dn=UnE)
UnS=0;                           %y=0 Neumann B.C (du/dn=UnS)
UnN=0;                           %y=L Neumann B.C (du/dn=UnN)

%%
%Initial Conditions
for i=1:nx
    for j=1:ny
        if ((0.5<=y(j))&&(y(j)<=1.5)&&(0.5<=x(i))&&(x(i)<=1.5))
            u(i,j)=2;
        else
            u(i,j)=0;
        end
    end
end

%%
% %B.C vector
% bc=zeros(nx-2,ny-2);
% 
% bc(1,:)=-UnW/dx;       %Neumann B.Cs
% bc(nx-2,:)=UnE/dx;     %Neumann B.Cs
% bc(:,1)=-UnS/dy;       %Neumann B.Cs
% bc(:,nx-2)=UnN/dy;     %Neumann B.Cs
% 
% %B.Cs at the corners: 
% bc=visx(x)*dt*bc+visy(y)*dt*bc;

%Calculating the coefficient matrix for the implicit scheme
Ex=sparse(2:nx-2,1:nx-3,1,nx-2,nx-2);
Ax=Ex+Ex'-2*speye(nx-2);
Ax(1,1)=-1; 
Ax(nx-2,nx-2)=-1;                 %Neumann B.Cs
Ey=sparse(2:ny-2,1:ny-3,1,ny-2,ny-2);
Ay=Ey+Ey'-2*speye(ny-2);
Ay(1,1)=-1; 
Ay(ny-2,ny-2)=-1;                 %Neumann B.Cs
A=kron(Ay/dy^2,speye(nx-2))+kron(speye(ny-2),Ax/dx^2);
D=speye((nx-2)*(ny-2))-visx(x)*dt*A-visy(y)*dt*A;
%%
%Calculating the field variable for each time step
i=2:nx-1;
j=2:ny-1;

for it=0:nt
    un=u;
    h=surf(x,y,u','EdgeColor','none');       %plotting the field variable
    shading interp
    axis ([0 2 0 2 0 2])
    title({['2-D Diffusion with alphax = ',visxstr];['2-D Diffusion with alphay = ',visystr];['time (\itt) = ',num2str(it*dt)]})
    xlabel('Spatial co-ordinate (x) \rightarrow')
    ylabel('{\leftarrow} Spatial co-ordinate (y)')
    zlabel('Transport property profile (u) \rightarrow')
    drawnow; 
    refreshdata(h)

    %Implicit method:
    % Source Term: Au(1-u/K)
    U=un;
    U(1,:)=[];
    U(end,:)=[];
    U(:,1)=[];
    U(:,end)=[];
    U=reshape(U+bc,[],1);
    U=D\U;
    U=reshape(U,nx-2,ny-2);
    u(2:nx-1,2:ny-1)=U;
%Neumann:
    u(1,:)=f(2,:)-UnW*dx;
    u(nx,:)=f(nx-1,:)+UnE*dx;
    u(:,1)=f(:,2)-UnS*dy;
    u(:,ny)=f(:,ny-1)+UnN*dy;
    pause(0.1);
end