Droplet Profile in Matlab- ODE stability

[Kudaros]

Hello,

I'm currently modeling the profile of a droplet (sessile drop, axisymmetric) in matlab. I've coded differential equations, applied the solver, and I get a reasonable result, except that it spirals continuously.

The ODE's in question are:
$$\frac{dx}{ds}=cos(\theta)$$
$$\frac{dz}{ds}=sin(\theta)$$
$$\frac{d\theta}{ds}=2b+cz-\frac{sin(\theta)}{x}.$$

ode23/45, etc produce similar results. I could limit my tspan so that the spiral doesn't continue, but I would like it to hit on a stable solution. I have yet to try a boundary value approach, but, while I am familiar with MATLAB in other areas, I have never used it to solve ODEs, so I would like to try this initial value approach first.

For solutions I use:
[s,x]=ode23(@(s,x)drpivp(s,x,p),sspan,x0);

where p contains two parameters and x0 contains initial angle theta, x, z values.
droplet ODE function:
function drpivp = drpivp(s,x,p);
%x(1)=theta
%x(2)=x
%x(3)=z
%r0 is curvature at apex
%c is capillarity constant
r0=p(1);
c=p(2);
drpivp=[2/p(1)+p(2)*x(3)-sin(x(1))/x(2); cos(x(1)); sin(x(1))];
%if you put an 'end' it reduces the required number of arguments - why?

Thanks!

 

[Kudaros]

Thinking about this further, I realized I initialized with $[\theta x z]=[1 1 1]$ just to avoid the NaN issue that rises when using the proper initialization of [0 0 0]. The set of equations is parametrized with respect to the arc length. At $s=0 ,x=z=\theta=0.$

To avoid this NaN issue, I need to change $\frac{d\theta}{ds}=b,$ as this is the condition at $s=0$. It isn't clear to me how to have the first ODE be something different for a particular value of s, though.