Attempting to solve 2 unknown functions using matrix method

My code below works, and its solving for temperature vs radial position along a cylindrical wall using nodes 1 through M with i being internal nodes within the walls boundaries. The energy generation, g, is a function of radius and that works perfectly.

The challenge is to also turn k into a function of position instead of being a constant, making two unknown functions to solve for at the same time. So solving for both temperature and k. If anyone could help me understand how to do that I'd appreciate it!

the k function is : k = 18.306 - (0.04246 * T) + (5.9733 e-5 * T^2) - (3.814 e-8 * T^3)

If anyone wants to copy & paste:

clear;
clc;
%%%%%INPUTS
Tf_in = 20;         % Fluid inside [C]
h_in = 100;         % HT Coeff inside [W/m^2K]
Tf_out = 100;       % Fluid outside [C]
h_out = 200;        % HT Coeff outside [W/m^2K]
k = 10;             % Thermal conductivity [W/mK]
r_in = .1;          % Internal diameter [m]
r_out = .2;         % Outer diameter [m]
L = r_out-r_in;     % Thickness of aluminum oxide cylinder wall [m]
a = 1e4;            % Thermal diffusivity [W/m^3]
b = 2e5;            % [W/m^4]
c = 5e7;            % [W/m^5]
M = 50;
Dr = L/(M-1);
r = r_in:Dr:r_out;
g = a + b*r+c*r.^2;
A = zeros(M,M);
B = zeros(M,1);
A(1,1) = (-h_in*2*pi*r_in)-k*(2*pi*r_in+Dr/2)/Dr;
A(1,2) = (k*2*pi*r_in+Dr/2)/Dr;
B(1) = -h_in*2*pi*r_in*Tf_in-g(1)*r_in*pi*r_in*Dr;
for i = 2:M-1
    A(i,i-1) = (k*(2*pi*r(i)-Dr/2))/Dr;
    A(i,i) = (-k*(2*pi*r(i)-Dr/2))/Dr-(k*(2*pi*r(i)+Dr/2))/Dr;
    A(i,i+1) = (k*(2*pi*r(i)+Dr/2))/Dr;
    g(i) = a+b*r(i)+c*r(i)^2;
    B(i) = -g(i)*2*pi*Dr*r(i);
end
A(M,M-1) = (k*(2*pi*r_out-Dr/2))/Dr;
A(M,M) = (-k*(2*pi*r_out-Dr/2))/Dr-h_out*2*pi*r_out;
B(M) = -g(M)*pi*r_out*Dr-h_out*2*pi*r_out*Tf_out;
T = A\B;
plot(r*1e2,T)
xlabel('Position (cm)')
ylabel('Temperature (C)')

Answer 1

[Torsten]

Now you have to solve a linear system A*T = B.

If k becomes temperature-dependent, your matrix A becomes dependent on T and you have to solve A(T)*T = B or equivalently A(T)*T-B = 0. Try "fsolve" to do this.

Do you try to solve a boundary value problem with your code ? Then it might be better to use "bvp4c" instead of discretizing the differential equation on your own.