matrix form initial condition for 2D pde

Hi I am new to MATLAB and really need your help. I have tried to write a code for a 2D parabolic PDE. first I discretized it using finite difference method (FTCS), then I wrote the following code (on college system). ICVF0 is the initial matrix (which is taken from rgb image) when I write the following code, it doesn't show any difference after specified time in fact the image shows the same image as the first one (initial condition):( . However it should be changed. I read MATLAB help several times but couldn't find any solution. I wonder if any one can help me on this? I appreciate any help in advance.

ICVFT=imread('D:\ICVFT.jpg');% these are true color images ICVF0=imread('D:\ICVF0.jpg'); SUV0=imread('D:\suv0.jpg'); SUVT=imread('D:\suvT.jpg'); %--------------------------------making them the same size----------------- [rowsA1,colsA1, numberOfColorChannelsA1] = size(ICVFT); [rowsA2, colsA2, numberOfColorChannelsA2] = size(ICVF0); [rowsB1, colsB1, numberOfColorChannelsB1] = size(SUV0); [rowsB2, colsB2 ,numberOfColorChannelsB2] = size(SUVT); % See if lateral sizes match. if rowsA2 ~= rowsA1 colsA1 ~= colsA2 % Size of B does not match A, so resize B to match A's size. ICVF0 = imresize(ICVF0, [rowsA1 colsA1]); end if rowsB2 ~= rowsA1 colsA1 ~= colsB2 SUVT = imresize(SUVT, [rowsA1 colsA1]); end if rowsB1 ~= rowsA1 colsA1 ~= colsB1 SUV0 = imresize(SUV0, [rowsA1 colsA1]); end %--------------------------------------------------------------------------

D=0.13; % diffusion coefficient in [mm^2 day^-1] alpha=2.3e-3; % alfa coefficient in [g^-1 ml day^-1] betha=1.9e-2; % betha coefficient in [day^-1]

L=159; % millimeters W=181; % millimeters [nx,ny]=size(ICVFT); dx=L/(nx-1); dy=W/(ny-1); total_t=input('please enter total time in days: '); %dt=input('please enter delta_t (time step size): '); dt=0.4; d1=D*dt/(dx^2); %diffusion number in x direction d2=D*dt/(dy^2); %diffusion number in y direction nt=1+(total_t/dt); % number of time steps Time=total_t; %-------------------------------------------------------------------------- %------------------------------- rho calculation---------------------------

A=alpha*SUV0-betha*ICVF0; B=ICVF0-(ICVF0.^2); rho0=A./B; AT=alpha*SUVT-betha*ICVFT; BT=ICVFT-ICVFT.^2; rhoT=(rho0+(AT./BT))/Time; %-------------------------------------------------------------------------- t=0; u=ICVF0; T=zeros(size(u)); for n=2:nt

 u=ICVF0;
 for i=2:nx-1
     for j=2:ny-1
         T(i,j)=d1*(u(i+1,j)+u(i-1,j))+d2*(u(i,j+1)+u(i,j-1))+(1-2*(d1+d2)+(rho0(i,j)+rhoT(i,j)*t)*dt*(1-u(i,j)))*u(i,j);
     end
 end
 T;
 t=t+dt;
  for i=2:nx-1
     for j=2:ny-1
        ICVF0(i,j)=T(i,j);
     end
 end
 ICVF0;
 end

figure,imtool(ICVF0);

[Jan]

The code is very hard to read. We do not have your image files, such that we cannot run your code. You do not ask a specific question. All we know is that you expect some changes, but do not observe this. So how could we create a useful answer?