I'm solving a set of 3 PDEs, 4 ODEs, and 1 AE in MatLab. The model describes the dynamic behavior of an argon/water vapor spherical bubble in water undergoing an oscillatory pressure field (P = P0+PA*sin(wt)). The first 2 PDEs describe the mass and energy conservation along the bubble radius [0,R], while the third PDE is the energy conservation in the liquid [R,+inf]:
The 4 ODEs are for the functions R(t), R'(t), P(t) and m(t), i.e. the time dependence of bubble radius, bubble wall speed, pressure in the bubble and total mass in the bubble. An additional AE is required to find the value of overdot m which is the mass flux at the interface. In order to solve the system, I've discretized the PDEs using a straightforward finite difference method with second-order accurate central differences for the spatial derivatives (properly modified at the boundaries to one-sided differences of the same accuracy order):
Therefore, I obtain a system of only ODEs and AEs that I solve defining a mass matrix and passing it to the ode15s solver. I am quite sure the code is properly written as it gives reasonable outputs for certain range of the parameters' value. However, for other ranges, the spatial derivatives at r=R of the liquid temperature PDE become very sharp and the code gives unrealistic outputs (or it just crashes before finishing the integration). In literature, it is advised to use an implicit formulation for that particular PDE. (For whoever works in a university there are good possibilities to have access to that paper at the following link PAPER. The following picture is an example of the spatial profile of the liquid temperature at different times obtained for certain parameters values.
If the acoustic pressure is increased above a certain threshold, the gradients of those profile at the bubble wall (r=R) should be much sharper (as the temperature goes up to several thousand Kelvins in few nanoseconds. That's when my code fails.
My question is whether it is possible to implement the implicit method within the tools I am using and, if yes, how to do it. I would like to keep using ode15s because the time dependence is also really steep and discretizing time as well would require writing a variable time step algorithm (which I am not familiar with).
Thanks in advance.