How can I plot this code?

Question

Hdez 2021-1-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/723208-how-can-i-plot-this-code

回答： Gaurav Garg 2021-1-25

clc
close all
clear all
%Lenght of the plate ----------------------------------------------
lx= input ('Enter the desire value for x length:');
ly= input ('Enter the desire value for y length:');
%alpha ------------------------------------------------------------
k=143;                                            %Thermal coductivity (W/m.C)
p=2.8*10^3;                                       %Density (kg/m^3)
c=795;                                            %Specific heat (J/kg.C)
alpha= k/(p*c);                                   %Thermal diffusivity (m^2/s)
%Grid Spacing
h= input('Enter the desire grid spacing:'); 
while rem (lx,h)~=0 || rem (ly,h)~=0
    fprintf ('Incorrect solution \n Please enter a value which results in an integer number:')
    h= input ('Reenter the desire grid spacing:');
end
%Grid Spacing along lx and ly ------------------------------------- 
nc= (lx/h)+1;                                     %Number of columns
nr= (ly/h)+1;                                     %Number of raws 
%Time step --------------------------------------------------------
dt= input ('Enter the desire time step:');       
while dt>= h^2/(4*alpha)
    fprintf('solution is unstable. Please return to dt and choose a smaller value for: %g\n', h^2/(4*alpha))
    dt= input ('Enter the desire time step:');
end
%Fourier Number (fo) ----------------------------------------------
fo= alpha*dt/(h^2);
%Creation of vectors ----------------------------------------------
T1= zeros (nr, nc);
%Boundary conditions ----------------------------------------------
T1(end,:)=0;
T1(1,:)=50;
for i= 1:nr
    T1(i,1)= 50-50/(nr-1)*(i-1)
    T1(i,end)= 50-50/(nr-1)*(i-1)
end
%For loops for temperatura: (Ts + Tn + Tw + Te - (4-1/fo)* Tp)) * fo
T2=T1;
for t=0:dt:inf
    T1=T2
    for j=2:nc-1
        for i=2:nr-1
            T2(i,j)= (T1(i-1,j)+ T1(i+1,j)+T1(i,j+1)+T1(i,j-1)-(4-1/fo)*T1(i,j))*fo;
        end
    end
    %check convergence 
    if max(max(T2-T1))<0.001
        break
    end   
end
disp(dt*t)