Solve equation with Crank Nicolson and Newton iterative method in Matlab

Question

I'm tring to solve the eqation something like this to get the transient temperature of each node:

where, is time; () is the input power; is the note ( = 1,2,.. .,n+1) (1 is the heat source and +1 is the constant temperature environment); = () is the temperature at while t; ℎ is the thermal resistance between and +1; ℎ is the thermal capacity of .

With the initial condition:

I try to discretize in time with Crank Nicolson method:

But I don't know how to do that with f(t) equation contains Ti-1, Ti and Ti+1 together, could anyone know how to solve this? Thanks a lot.

n = 4; 
dt = 0.01; 
timesteps = 1000; 
R = [6.9364, 8.3004, 2.8952, 12.0162]; 
C = [1, 2, 1, 2]; 
P_steady_state = 1.24; 
T_steady_state = P_steady_state * sum(R); 

T = 25*ones(n, timesteps);
T(:, 1) = T_steady_state;

% Crank-Nicolson method
for t = 2:timesteps-1
    for j = 1:n-1
        if j == 1
           f = P_steady_state + T(j+1, t) - T(j, t) / (C(j) * R(j));
        else
           f = (T(j-1, t) - T(j, t))/(C(j) * R(j-1)) + (T(j+1, t) - T(j, t))/(C(j) * R(j));
        end
        if j == 1
           f_t = P_steady_state + T(j+1, t+1) - T(j, t+1) / (C(j) * R(j));
        else
           f_t = (T(j-1, t+1) - T(j, t+1))/(C(j) * R(j-1)) + (T(j+1, t+1) - T(j, t+1))/(C(j) * R(j));
        end
        T(j, t+1) = T(j, t) + 0.5 * dt * f_t + 0.5 * dt * f;
    end
end

I have tried something like this, but it looks wrong.