% 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');