Inequality boundary condition (w<inf) and passing condition between solvers (w1 == w2)

Hello, I am failing in finding solution for a following problem:

I have a system of three consecutive differential equations, where each solves deformation "w" (in z direction) depending on radius "r" of an element which has 1-3 layers (piezoelectric sandwich) where each segment (1,2,3 layers) have its own similar DE. I need somehow tell the solver (currently ode45) two things about the conditions:

1) that at r=0, the calculated variable w<inf (because of the division by zero) - it only accepts numbers

2) each interval of the deformed element is connected by inner boundaries, i.e. w1 = w2 (condition for one DE is the solution of other DE), or Dw1 = Dw2

I want to do it purely numerically, so I don't want to express the conditions from equation, because I have an analytical solution and FEM solutions and need to compare them with the numerical solution. This also means that I know the system, I even know the solution, I just need to tell matlab to recalculate it for me numerically.

I tried other solvers, I tried checking BVPs. The issue is, that the only boundary condition is that at the ends of the element the deformation is zero, all of the rest is calculated and the origin of the deformation is from within the element itself (piezoelectric effect).

Would it be better to use simulink for this?

*Based on request, including the equations in question (although the question is general):

And conditions on the boundaries:

and between segments (only one set of them):

I just need to somehow put the inner and inequality conditions instead of the zeroes:

[Sol3] = ode45(f3,[0 R_s],[0, 0, 0, 0, 0]);
[Sol2] = ode45(f2,[R_s R_a],[0, 0, 0, 0, 0]);
[Sol1] = ode45(f1,[R_a R_m],[0, 0, 0, 0, 0]);

[Torsten]

It would help if you included the equations you want to solve together with the boundary conditions.

[Vojtech Lindauer]

Included it, although I don't think it is necessary, the question is general.

[Torsten]

You need five boundary conditions, but you only specified three. The conditions of finiteness at r = 0 do not count.

From heat conduction problems in a cylinder or in a sphere, the usual condition at r = 0 is dT/dr = 0 for the temperature T. Maybe it's similar here ?

[Vojtech Lindauer]

I need 15 (3 DEs * 5 variables) and I have 15, including at r=0. As I said I have analytical solution. The calculation is quite... large... And I only have a particular issue. So why clutter the question with all of it. I would have to include all the geometries, all the coefficient etc. All I need is to connect the solvers so that they work together on the entire interval <0,R_m>.

[Vojtech Lindauer]

Well yeah, that make sense given it is the maximum, for the Dw to be zero. I could use this workaround for the numerical solution, in the analytical solution these two conditions at r=0 are quite important (they rid me of 3 terrible coefficients :D). That still does not solve the issue of connecting the each DE at each interval though.

[Torsten]

The coupling at the interior points is done by exploiting the feature of bvp4c to solve Multipoint Boundary Value problems:

https://uk.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

As said: Every numerical solver needs 5 explicit exterior boundary conditions in sum. A condition that only says that a certain function value or derivative is finite is useless for a numerical computation. In heat transfer problems, the condition that T is finite at r = 0 in cylindrical or spherical coordinates is replaced by dT/dr = 0 at r = 0, e.g. This will numerically exclude log(r) or 1/r terms that appear in the general analytical solution.

[Vojtech Lindauer]

Yes, that is exactly what happens, it gets rid of terms with 1/r. The page you linked seems like something I was desperate to find, but failed to do so. Will try.