Trouble with my for loop and graph

The following code is not plotting the temperature distriutuon as desired. The graph should be much smoother. I ran this by my teacher. He said all of my coeficients are good, but was no help in fixing the code. All the math is good, worked out and approved. something is just weird with my code. Hoping someone can see something obvious that I may be missing. Edit: I added an attachment of the i=2:n-2 loop, as well as the governing equations/PDE. The boundary conditions are just the Ri and Rf values with the initial conditions of Ti= 20 degrees C and Tf is 790 C

clear
clc
close all

%Define Constants
rho = 150;
cp = 2300;
k = 0.30;
Ri = 0.0716;
Ro = 0.0970;
Ti = 20; 
Tf = 790;
alpha = k/(rho*cp);
q_ppp = -100000;

%Position and Time
N=51;
r=linspace(Ri, Ro, N);
dr=r(2)-r(1);
ti=0;
tf=60;
dt=0.01;
t=ti:dt:tf;
M=length(t);

T=zeros(N,M); %(Row,Column)
T(:,1) = Ti; %Ti
T(N,:)= Tf; %Tf

% % 5 Spatial Discretization
A = zeros(N-1,N-1);
d = zeros(N-1,1);
b1 = -(2*alpha/(dr^2)); % for i=1, T2
c1 = 2*alpha/(dr^2); %i=1, T1
bn = -2*alpha/dr^2; %bn-1, %Tn-1
an = ((alpha/(dr^2))*(1-(dr/(r(N-1))))); %an-1, Tn-2
A(1,1) = b1;
A(1,2) = c1;
A(N-1,N-1) = bn;
A(N-1,N-2) = an;

for j = 2:N-2
    A(j,j) = (-2*alpha/(dr^2)); %b, Ti
    A(j,j-1) = (alpha/dr^2)-(alpha/(r(j)*dr)); %a, %ti-1
    A(j,j+1) = (alpha/dr^2)+(alpha/(r(j)*dr)); %c, Ti+1
end
d(1) =  q_ppp/(rho*cp);
d(N-1) = q_ppp/(rho*cp);

for j = 2:N-2
    d(j) = q_ppp/(rho*cp);
end
% 6 Time Discretization
 
I = eye(N-1);
B = I - (dt/2)*A;
C = I + (dt/2)*A;
Binv = inv(B);

for n=1:M-1
   T(1:N-1,n+1) = Binv*C*T(1:N-1,n) + Binv*dt*d;
end


disp('Plot')

plot(r,T(:,3000),'b')
title('Temperature Distribution')
xlabel('X')
ylabel('y')

hold on
grid on
grid minor 
box on 

[Torsten]

Could you state that PDE including initial and boundary conditions you are trying to solve ?

One obvious error in your discretization is that you don't reference T(N,:) in your matrix A. This should at least be the case in the equation where you update T(N-1,:), but it isn't.

[Kyle Langford]

@Torsten I added some attachments of the problem statement.

[Torsten]

Could you reformat the code ?

[Torsten]

Here is a code you can compare your results with.

As I said: You never reference the outward temperature Tf in your code so as to initialize heat conduction from the boundary.

%Define Constants
rho = 150;
cp = 2300;
k = 0.30;
Ri = 0.0716;
Ro = 0.0970;
Ti = 20; 
Tf = 790;
alpha = k/(rho*cp);
q_ppp = -100000;

x = linspace(Ri,Ro,501);
t = linspace(0,60,100);
m = 2;
sol = pdepe(m,@(x,t,u,dudx)heatsphere(x,t,u,dudx,rho,cp,k,q_ppp),@(x)heatic(x,Ro,Ti,Tf),@(xl,ul,xr,ur,t)heatbc(xl,ul,xr,ur,t,Tf),x,t);
u = sol(:,:,1);

plot(x,[u(10,:);u(30,:);u(50,:);u(70,:);u(90,:);u(100,:)])
grid on

function [c,f,s] = heatsphere(x,t,u,dudx,rho,cp,k,q_ppp)
  c = rho*cp;
  f = k*dudx;
  s = q_ppp;
end

function u0 = heatic(x,Ro,Ti,Tf)
  if x==Ro
    u0 = Tf;
  else
    u0 = Ti;
  end
end

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t,Tf)
  pl = 0; 
  ql = 1; 
  pr = ur-Tf;
  qr = 0; 
end

[Kyle Langford]

@Torsten Yes, sorry about that. I didn't realize it went back to text format.

I am not sure where to reference the T matrix in my A matrix. I am supposed to be finding dT/dt=A*T+d, but we don't know T besides Ti and Tf.

After developing the time discretization Solving for T gives B^-1*C*T+B^-1*dt*d

I am just confused on how to apply this to the code.

[Kyle Langford]

@Torsten That's beautiful. I willt try to manipulate my code to get similar results. that is really clean. I am still pretty new to coding, so I don't know all the best ways to manipulate code.

Thank you for the help!