Heat flow in 1d

Yes, you are right, because θ is an array which is mutable (we can modify its contents) so when you run the code for 2nd time, it would use the θ values from 1st iteration (which are not zeros anymore) and so on…
so probably u just need to update the whole array once in one loop and then you can plot the complete θ values just once, which would look like exponential decreasing scatter plot.
Code is as follows →

using Plots

begin
Δg=0.000000001 ;m=90
b=-2-(m^2)*(Δg)^2
end

θ=zeros(100)

begin
θ[1]=1;θ[100]=0
end

p=plot();

for i in 2:99
    θ[i] = round((-1/b)*(θ[i-1] + θ[i+1]),digits=9);
    print(θ[i]);
end

plot(θ[:], seriestype=:scatter)