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