% Stern-wave toy problem
clear
close all
ep = 0.5; % Try different values of eps from 0.2 to 1
f0 = 1e-3; % Set the domain from [f0, fmax]
fmax = 50;
% Leading-order solution
q0 = @(f) (f./(f+1)).^(1/2);
% Define the differential equation
dq = @(f, q) (-1i/2)*(q.^2 - q0(f).^2)./(ep*q^2.*q0(f));
[phi, Q] = ode45(dq, [f0, fmax], q0(f0));
qsol = Q(:,1);
figure(1);
plot(phi, real(qsol), 'b');
hold on
plot(phi, q0(phi), 'r');
hold off
xlabel('\phi'); ylabel('q(\phi)');
title('Solution of the model stern-wave problem');