Piecewise ODE second order

Hello,

I am trying to solve the second order ODE as seen below:

This corresponds to the conductive heat transfer and the Joule heating of the profile seen 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 Tk, Tc and Tf, respectively. (The subscripts stand for hot-beam, cold-beam and flexure)

So I want to solve the ODE three times

The six boundary conditions that must be taken into account can be seen below:

The solution must be a continuous function and should look something like this:

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!

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 on the picture boundary conditions (temperatures)?