I'm not sure how long it took you, but your code as it's written only took 1-2s, which wouldn't make a good problem to parallelize. I'm assuming n is much larger and that the code takes much longer to run.
With that said, I'm not sure this can be parallelized. To begin with, I'll show you how I would resolve the classification of a (and there are others, such as b and c). In doing so, you'll see how the parfor has a loop dependence on p0.
I've also added a few suggestions to cleanup/vectorize your code (listed as "% MW") -- there's more that could be done, but this gives you the idea.
function ismailov
n=91; h=1/(n-1); nt=1080; tau=0.001;
% MW: Preallocate x & y
x = nan(n,1);
y = nan(n,1);
for i=1:n
x(i)=(i-1)*h; y(i)=(i-1)*h;
end
% MW: No need for double for-loop
% for j=1:n
% for i=1:n
% p0(i,j)=1;
% end
% end
p0 = ones(n);
t = nan(nt,1);
pt = nan(n);
for jt=1:nt
t(jt)=jt*tau;
for j=1:n % Anlitik yechimni hisoblash
for i=1:n
pt(i,j)=exp(-2*x(i)*y(j)*t(jt));
end
end
end
tic
parfor jt=1:nt % Vaqt bo`yicha tsikl}
p0 = unit_of_work(n, h, tau, jt, t, x, y, p0);
end % jt
toc
% MW: No need for double for-loop
% for j=1:n
% px(j)=p0(21,j);
% pxt(j)=pt(21,j);
% end
px = p0(21,:);
pxt = pt(21,:);
x=0:h:1;
y=0:h:1;
figure ('Position',[600 160 700 600]);
subplot(3,2,1);
meshc(x,y,pt) %Uch o'lchovli 3D grafika
title('АНИК ЕЧИМ 3D ГРАФИГИ')
subplot(3,2,2);
meshc(x,y,p0) %Uch o'lchovli 3D grafika
title('CОНЛИ ЕЧИМ 3D ГРАФИГИ')
subplot(3,2,3);
contour(x,y,pt,'ShowText','on','LineWidth',2); %kontur chizish
title('АНИК ЕЧИМ КОНТУР ГРАФИГИ')
subplot(3,2,4);
contour(x,y,p0,'ShowText','on','LineWidth',2); %kontur chizish
title('CОНЛИ ЕЧИМ КОНТУР ГРАФИГИ')
subplot(3,2,5);
plot(x,px,x,pxt); %Ikki o'lchovli grafika
title('КЕСИМДА АНИК ВА CОНЛИ ЕЧИМ ЎЗГАРИШИ X БЎЙИЧА')
end
function p0 = unit_of_work(n, h, tau, jt, t, x, y, p0)
% MW: Preallocate a, b, c, d, aa1, bb1, p
a = ones(1, n-1);
c = ones(1, n-1);
b = rand(1, n-1);
d = rand(1, n-1);
aa1 = rand(1, n-1);
bb1 = rand(1, n-1);
p = rand(n);
for j=2:n-1 % Tenglamalar koeffitsientini hisoblash
for i=2:n-1
% a(i)=1; c(i)=1; b(i)=2+2*h*h/tau;
b(i)=2+2*h*h/tau;
ff=h*h*(-2*(2*t(jt)*t(jt)*(x(i)*x(i)+y(j)*y(j))+x(i)*y(j))*exp(-2*x(i)*y(j)*t(jt)));
d(i)=2*h*h/tau*p0(i,j)+(p0(i,j-1)-2*p0(i,j)+p0(i,j+1))+ff;
end
aa1(1)=0; bb1(1)=1;
for i=2:n-1
aa1(i)=c(i)/(b(i)-aa1(i-1)*a(i));
bb1(i)=(d(i)+bb1(i-1)*a(i))/(b(i)-aa1(i-1)*a(i));
end
p(n,j)=exp(-2*y(j)*t(jt));
for i=n-1:-1:1
p(i,j)=aa1(i)*p(i+1,j)+bb1(i);
end
end %j
for i=1:n % Chegaraviy qiymatlarni hisoblash
% p(i,n)=(4.0*p(i,n-1)-p(i,n-2)+h)/3; % p(i,1)=(4.0*p(i,2)-p(i,3)+2*h)/3;
p(i,n)=exp(-2*x(i)*t(jt)); p(i,1)=1;
end
% MW: No need for double for-loop
% for j=1:n % k+1 qatlam uchun boshlangich qiymatni hisoblash
% for i=1:n
% p0(i,j)=p(i,j);
% end
% end
p0 = p;
for i=2:n-1 % Tenglamalar koeffitsientini hisoblash
for j=2:n-1
a(j)=1; c(j)=1; b(j)=2+2*h*h/tau;
ff=h*h*(-2*(2*t(jt)*t(jt)*(x(i)*x(i)+y(j)*y(j))+x(i)*y(j))*exp(-2*x(i)*y(j)*t(jt)));
d(j)=2*h*h/tau*p0(i,j)+(p0(i-1,j)-2*p0(i,j)+p0(i+1,j))+ff;
end
aa1(1)=0; bb1(1)=1;
for j=2:n-1
aa1(j)=c(j)/(b(j)-aa1(j-1)*a(j));
bb1(j)=(d(j)+bb1(j-1)*a(j))/(b(j)-aa1(j-1)*a(j));
end
p(i,n)=exp(-2*x(i)*t(jt));
for j=n-1:-1:1
p(i,j)=aa1(j)*p(i,j+1)+bb1(j);
end
end %i
for j=1:n % Chegaraviy qiymatlarni hisoblash
%p(n,j)=(4.0*p(n-1,j)-p(n-2,j)+2*h)/3; %p(1,j)=(4.0*p(2,j)-p(3,j)+2*h)/3;
p(n,j)=exp(-2*y(j)*t(jt)); p(1,j)=1;
end
% MW: No need for double for-loop
% for i=1:n % keyingi qatlam uchun boshlangich qiymatni hisoblash
% for j=1:n
% p0(i,j)=p(i,j);
% end
% end
p0 = p;
end
Let's look closer at the parfor loop
parfor jt=1:nt % Vaqt bo`yicha tsikl}
p0 = unit_of_work(n, h, tau, jt, t, x, y, p0);
end % jt
I've refactored the for-loop to call unit_of_work. The only variable modified in your for-loop is p0. But notice that p0 is not dependent on the loop variable jt. For instance
parfor jt=1:nt
p0(jt,:) = unit_of_work(n, h, tau, jt, t, x, y);
end % jt
In this case, we would be building up each row independent of any other row, which then allows us to assign each row to a worker. Instead, the entire matrix p0 is dependent on the previous state of the matrix -- what I'd call a cascading algorithm -- and therefore can't be parallelized.
