Plotting Only 4 Time steps

Hello,

I really failed to make my plot draw only the results of 4 time spots out of the whole time steps, I tried to use hold on but it drew every single iteration. I need just to show 4 screen shots of the evoloving temperature disutribution (begining, middle, middle 2, final), can anyone help me pls.

Thank you

[the cyclist]

Here is the code, for those who do not care to download it ...

%Seconf Deravitive solution

set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',17)
clear *
L=0.025; % copper wire length
N = 100; % you decide, # of discretization points
% the more you increase the N number, the more it is close to the BC
NN = N + 2; % 2 represent the ghost cells.
dx = L/N;
x = linspace (-dx/2,L+dx/2,NN) % define temp. solution dx/2 3dx/2 etc..

T = 293 + sin(2*pi*x/L);


in = 2:1:(NN-1); % Internal cell
rhs = in+1 % Right hand side
lhs = in-1 %Left hand side

d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2); %1st derivatve of the heat equation shown in the slide page 37
Ta = -sin(2*pi*x/L)*(2*pi/L)^2; %analytical

T = 293 +sin(2*pi*x/L);
%Ta = cos(2*pi*x/L)*(2*pi/L); % analytical derivative (maybe)
figure(2),clf
plot(x(in),d2Tdx2,'o-' ), hold on% we have put n to make the vectors match by excluding the ghost cells.
plot(x(in),Ta(in),'r-')

legend('Numerical d2Tdx2','Analytical dTdx')

%%%%%%%%%%%%%%%Start the HW%%%%%%%%%%
%solve dTdt = alpha*d2Tdx2
%using Euler method and central differnce
% for space derivative

alpha = 10.1; %m2/s
T = 300*ones(1,102)% Initalize
Tmin = 310; Tmax=300;
T(1) = 2*Tmin - T(2) %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);

dt = 0.000000001; %timestep in second

simutime = 1; %seconds
simusteps = round(simutime/dt);

for(t=1:simusteps)

d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2);
dT = alpha*dt*d2Tdx2; %get dT in each cell
T(in) = T(in) + dT;
%set BC
%If you want to isolate the left side
% you can add T(1) = T(2)
% and comment the line below.
T(1) = 2*Tmin - T(2) %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);

if(mod(t,10))
    figure(3);
    plot(x(in),T(in),'r-','LineWidth',1), hold on
    h=xlabel('Length (m)');
    h=ylabel('Temperature (K)');
    title('Part A');
end
end 

Answer 1

[the cyclist]

I expect that this if statement is not doing what you expect ...

if(mod(t,10))

The lines following that will be executed whenever mod(t,10) is non-zero, which is most iterations.

[Fares]

Thank you very much,

I need to make my graph capture 4 iterations, just like this figure. I really dont know how to make it capture few iterations while it calculates...

[the cyclist]

I made a few changes to your code. The most important is the change to that mod() function, so now it only plots every 5,000 iterations.

Note that your code runs for a billion iterations.

set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',17)
clear *
L=0.025; % copper wire length
N = 100; % you decide, # of discretization points
% the more you increase the N number, the more it is close to the BC
NN = N + 2; % 2 represent the ghost cells.
dx = L/N;
x = linspace (-dx/2,L+dx/2,NN); % define temp. solution dx/2 3dx/2 etc..

T = 293 + sin(2*pi*x/L);


in = 2:1:(NN-1); % Internal cell
rhs = in+1; % Right hand side
lhs = in-1; %Left hand side

d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2); %1st derivatve of the heat equation shown in the slide page 37
Ta = -sin(2*pi*x/L)*(2*pi/L)^2; %analytical

T = 293 +sin(2*pi*x/L);
%Ta = cos(2*pi*x/L)*(2*pi/L); % analytical derivative (maybe)
figure(2),clf
plot(x(in),d2Tdx2,'o-' ), hold on% we have put n to make the vectors match by excluding the ghost cells.
plot(x(in),Ta(in),'r-')

legend('Numerical d2Tdx2','Analytical dTdx')

%%%%%%%%%%%%%%%Start the HW%%%%%%%%%%
%solve dTdt = alpha*d2Tdx2
%using Euler method and central differnce
% for space derivative

alpha = 10.1; %m2/s
T = 300*ones(1,102);% Initalize
Tmin = 310; Tmax=300;
T(1) = 2*Tmin - T(2); %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);

dt = 1.e-9; %timestep in second

simutime = 1; %seconds
simusteps = round(simutime/dt);

for t=1:simusteps
    
    d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2);
    dT = alpha*dt*d2Tdx2; %get dT in each cell
    T(in) = T(in) + dT;
    %set BC
    %If you want to isolate the left side
    % you can add T(1) = T(2)
    % and comment the line below.
    T(1) = 2*Tmin - T(2); %0.5*(T(1)+T(2)) = Tmin
    T(NN) = 2*Tmax - T(NN-1);
    
    figure(3)
    hold on
    if mod(t,5e3)==0
        plot(x(in),T(in),'-r','LineWidth',1)
        xlabel('Length (m)');
        ylabel('Temperature (K)');
        title('Part A');
    end

end

[Fares]

That helps a lot.

Thank you very much for your help, it is most appreciated.