ode45 errors in integration

I'm trying to solve a second-order differential equation (Young-Laplace) with ode45() using shooting method. The solution should be a function u(z), evaluated in the interval [zi zf], that respects the initial and final boundary conditions u(zi) = bci and u(zf) = bcf.

To find the correct solution, the ODE needs to be split into two 1st order equations that ode45() can integrate along. By providing u(zi) = bci and trying different values of u'(zi), I can converge on the initial slope that will give the correct u(zf). This summarizes the shooting method.

The issue is that for some reason the value of u(zf) is not a smooth monotonous function of u'(zi) when integrating with ode45(). See figure below:

I believe the solutions of the differential equation should provide a smooth function. Which means ode45() is having some issue.

In the case depicted the target final boundary u(zf) is marked by the horizontal black-line. And as can be seen from the inset the non-smoothness is preventing my shooting method from converging on a solution, because a glitchy singularity happens to be located around the region.

Is there an aleternative to ode45()? Or am I using it wrongly?

The code illustrating the issue and used to create the figure above is shown below:

clear all; 
clc;
bci = 3.84915307352372e-05 
bcf = 0.0005127;
H = -2791.52777777778;
zi = 0;
zf = 0.00122121172192511;

u = @(z,u) [u(2);  (1 + u(2)^2)/u(1) + H* (1 + u(2)^2)^(3/2)]; % Young-Laplace equation

slopes = 2:.01:15; % Values of slope to be tested.
for i = 1:length(slopes)
    [z, r] = ode45(u, [zi zf], [bci slopes(i)]);
    lastValue(i) = r(end,1);
end

figure(1); % u(zf) as function of u'(zi).
clf;
hold on;
plot(slopes, lastValue, '.-');
yline(bcf);
hold off;
xlabel('Slope at z_i  - u''(z_i)');
ylabel('Last value of integration - u(z_f)');

figure(2);
clf;
hold on;
plot(z.*10^6, r(:,1).*10^6, '.-');
plot([zf zf ].*10^6, [0 bcf].*10^6, '.-', 'Color', 'r'); 
plot([zi zi].*10^6, [0 bci].*10^6, '.-', 'Color', 'r');
hold off;
xlim([-20 1300]);
ylabel('u(z)');
xlabel('z');
title('Solution for last slope tested');



%% Plot inset
figure(3);
clf;
hold on;
plot(slopes, lastValue, '.-', 'LineWidth', 3,'MarkerSize',22);
yline(bcf);
hold off;
xlabel('Slope at zi  - u''(z_i)');
ylabel('Last value of integration - u(z_f)');
xlim([7.79 8.45]);
ylim([512.029 513.687]*10^-6);

I do not show the code for the shooting method, as I believe the issue is well explained with this alone.

Let me know if more information is needed.

Answer 1

[Torsten]

Try setting RelTol and AbsTol to smaller values than the defaults:

options = odeset('RelTol',1e-6,'AbsTol',1e-8);
[z,r] = ode45(...,options)

[Arthur Vieira]

That seems to have done the trick. Thanks!