MATLAB 偏微分方程式　数値シミュレーション

MATLABを用いた、2次元偏微分方程式(反応拡散方程式)の濃度変化のシミュレーションについて質問します。

clear all;

%*各種変数の設定(時間刻み幅、格子間隔など)。

N = input('格子点の数を入力：');
L = 1.; % 系の辺
h = L/(N-1); % 格子間隔
time = input('時間刻み幅を入力：');
D = 10^(-4); % 拡散係数
coeff = D * time / h^2;
alpha = 3.0*10^(-2); % the late of release
beta = 1; % the stlength of cell responce
tau = 64*60; % the characteristic degradation time of chemoattractant

%* 初期条件と境界条件を設定

nn = randn(N,N);
nn_new = zeros(N, N);
vvx = zeros(N, N);
vvy = zeros(N, N);
vvx_new = zeros(N, N);
vvy_new = zeros(N, N);
cc = zeros(N, N);
cc_new = zeros(N, N);

xplot = (0:N-1) * h - L/2;
yplot = (0:N-1) * h - L/2;
iplot = 1;
nstep = input('時間刻みの数を入力：');
nplots = 50;
plot_step = nstep/nplots;

clf;   

%差分法                                                             
    
for istep = 1:nstep
        for i = 2:N-1
            for j = 2:N-1
            cc_new(i,j) = cc(i,j) + ...
                    coeff*(cc(i+1,j) + cc(i-1,j) + cc(i,j+1) + cc(i,j-1) -4*cc(i,j)) + ...
                time*(alpha * nn(i,j) - (1/tau) * cc(i,j));
        end
    end
    
    for i = 2:N-1
        for j = 2:N-1
            vvx_new(i,j) = vvx(i,j) - (time/(2*h)) * ((vvx(i+1,j) - vvx(i-1,j)) * vvx(i,j) + ...
                (vvy(i+1,j) - vvy(i-1,j)) * vvy(i,j)) + ...
                (time/(2*h)) * beta * (cc(i+1,j) - cc(i-1,j));
            
            vvy_new(i,j) = vvy(i,j) - (time/(2*h)) * ((vvy(i,j+1) - vvy(i,j-1)) * vvx(i,j) + ...
                (vvy(i,j+1) - vvy(i,j-1))*vvy(i,j)) + ...
                (time/(2*h)) * beta * (cc(i,j+1) - cc(i,j-1));
        end
    end
    
    for i = 2:N-1
        for j = 2:N-1
            nn_new(i,j) = nn(i,j) + ...
                (time/(2*h)) * (nn(i+1,j) * vvx(i+1,j) - nn(i-1,j) * vvx(i-1,j) + ...
                nn(i,j+1) * vvy(i,j+1) - nn(i,j-1) * vvy(i,j-1));
        end
    end
    
    nn = nn_new;
    
    vvx = vvx_new;
    
    vvy = vvy_new;
    
    cc = cc_new;
    
    nn
    
    if(rem(istep,plot_step) < 1)
        nnplot = nn;
        iplot = iplot + 1;
        figure(iplot);
        mesh(xplot,yplot,nnplot);
        xlabel('x');     ylabel('y'); zlabel('N(x,y)');
        grid on;
        hold on;
    end
end