How can I solve a piecewise ODE

I am trying to solve the second order ODE as seen above, which corresponds to the conductive heat transfer and the Joule heating of the profile seen below. As seen in the figure below, the profile has varying thickness causing the parameters with the subscript 'i' to change for each of the intervals.

The temperature in the beams with length L1, L2 and L3 will now be denoted as , and , respectively. (The subscripts stand for hot-beam, cold-beam and flexure)

The six boundary conditions that must be taken into account are stated below and are indicated in the figure above:

(1)

(2)

(3)

(4)

(5)

(6)

The solution must be a continuous function and should look something like the graph seen below.

So far I have tried to solve this problem using the function dsolve and using each interval as a seperate ODE. This did however not result in the desired solution, as the boundary conditions are not met and the resulting graph is not a continuous function. The lines of code that I used have been added below.

I would like to know what I should do different or which function to use to solve this problem. Thanks in advance!

Dion

clear all; close all; clc;
%% Constants
    % Dimensions
    
h=10e-6;                        % height
b=[20e-6 60e-6 20e-6];          % width
L=[200e-6 140e-6 40e-6];        % length
P=2*h+2*b;                      % perimeter
A=h*b;                          % area

    % Thermal properties
    
U=15000;        % conduction heat transfer coefficient
k_p=50;         % thermal concuctivity
T_inf=293.15;   % ambient temperature
V=[5 5 5];      % applied voltage
r=22e-6;        % resistivity
R=r*L./A;       % electrical resistance



%% Differential Equations
syms T_h(x) T_c(x) T_f(x)
DT_h = diff(T_h);
DT_c = diff(T_c);
DT_f = diff(T_f);

ode1 = diff(T_h,x,2) == ((U*P(1))/(A(1)k_p))(T_h-T_inf)-(V(1)^2/(A(1)*k_p*R(1)*L(1)));
ode2 = diff(T_c,x,2) == ((U*P(2))/(A(2)k_p))(T_c-T_inf)-(V(2)^2/(A(2)*k_p*R(2)*L(2)));
ode3 = diff(T_f,x,2) == ((U*P(3))/(A(3)k_p))(T_f-T_inf)-(V(3)^2/(A(3)*k_p*R(3)*L(3)));

odes=[ode1 ode2 ode3];
%% Conditions
cond1 = T_h(0) == 293.15;
cond2 = T_f(L(3)) == 293.15;
cond3 = T_h(L(1)) == T_c(0);
cond4 = T_c(L(2)) == T_f(0);
cond5 = A(1)*DT_h(L(1)) == A(2)*DT_c(0);
cond6 = A(2)*DT_c(L(2)) == A(3)*DT_f(0);

conds = [cond1 cond2 cond3 cond4 cond5 cond6];
%% Solve
[T_h(x), T_c(x), T_f(x)]=dsolve(odes,conds);
%% Plotting
x1=linspace(0,L(1),1000);
x2=linspace(0,L(2),1000);
x3=linspace(0,L(3),1000);

xtot=[x1 L(1)+x2 L(2)+L(1)+x3];

T1=subs(T_h,x,x1);
T2=subs(T_c,x,x2);
T3=subs(T_f,x,x3);

Ttot=[T1 T2 T3];

figure()
plot(xtot, Ttot)
grid on

[darova]

Can you show points ( ) and temperatures (boundary conditions) on the picture?

[Dion Stam]

I'm sorry for the late reply. I numbered the boundary conditions and showed their respective position in the top figure. I hope this clears things up.

[darova]

WHere are points ?

[Dion Stam]

h, c and f are referring to the hot, cold and flexure part of the structure, respectively. This is illustrated in the graph. The hot arm is the interval (0,L1), the cold arm is between (L1,L2) and the flexure is between (L2, L3).

[darova]

Can you explain (5) and (6) conditions. Are those related to the whole region or some points?

[Dion Stam]

Conditions (5) and (6) imply that the heat transfer rate must be constant throughout the entire structure, but are applied at L1 and L2 since the area changes at these points.

[darova]

I don't quite understand these conditions

They mean that form of curves are the same. But ,

[Dion Stam]

The boundary conditions are obtained by utilising the continuity of both temperature T, and the rate of heat conduction across the joint points of the hot arm, cold arm and flexure. I used these conditions based on a previous study on this subject. If you are interested I added the hyperlink to the study below (page 8 shows the boundary conditions), but I understand if this is a little too much.

https://www.mdpi.com/2072-666X/9/3/108/htm

[darova]

Sorry im lost. I don't know how to proceed

Answer 1

[Emilien Flayac]

If this can still be useful, you could try to put your problem into the following form:

where would be evaluated on grids and Fwould contain your dynamical constraint as well as the boundary conditions (1)-(6).

Then, you could use any method to find the roots. Since you already used the Symbolic toolbox, have a look at vpasolve.