Метод явных сеток для уравнения диффузии

% Метод сеток для уравнения диффузии
c=[];x=[];
n=9;                                         % желаемое число кривых при разных t
q=2.^n;                                    %  число шагов по времени
h=1./20;                                   % шаг по x
d=.05;                                      % коэффициент диффузии 
t=.01;                                       % шаг по времени
a=d.*t./h.^2;                            % параметр
v=0;                                       % скорость конвективного переноса
w=v.*t./h;                                 % параметр
% Начальные условия (выберите, измените)
c(1,:)=[0 0 0 0 0 0 0 .2 .2 .2 .2 .2 .2 .2 0 0 0 0 0 0 0];
% c(1,:)=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1 .9 .8 .7 .6 .5 .4 .3 .2 .1 0];
for j=1:21
x(j)=(j-1)./20;
end
for i=1:q                                 % время
for j=2:20                               % положение по х
% Варианты граничных условий - выберите
c(i,1)=0;c(i,21)=0;
% c(i,1)=c(i,2)./(1+v.*h./d);c(i,21)=c(i,20)./(1-v.*h./d);
% c(i,1)=0;c(i,21)=c(i,20)./(1-v.*h./d);
% c(i,21)=0;c(i,1)=c(i,2)./(1+v.*h./d);
% Введение внешнего источника в точку j=1
if j==11
m(i,j)=0;                                % введите 1, если есть
else m(i,j)=0;end
% Собственно процедура расчета
c(i+1,j)=-w.*(c(i,j+1)-c(i,j))+c(i,j)+a.*(c(i,j+1)-2.*c(i,j)+c(i,j-1))+m(i,j).*t;
end;end
% Графический вывод результатов
figure(1)
plot(x,c(1,:));hold on;axis([0 1 0 .25]);pause(2);
for k=1:n
p=2.^k;t1=t.*p;
plot(x,c(p,:),'r');axis([0 1 0 .25]);
pause(1);
% hold off                      % удалите  "%", если удаляем предыд. кривые
hold on                          % вставьте "%", если удаляем предыд. кривые
end
xlabel('position')
ylabel('concentration')
hold off

[X,T] = meshgrid(x,t);
figure
surf(X,T,c,'facealpha',0.7)
grid on
shading interp
xlabel('position')
ylabel('time')
zlabel('concentration')