1. % Stern-wave toy problem
  2.  
  3. clear
  4. close all
  5.  
  6. ep = 0.5;       % Try different values of eps from 0.2 to 1
  7. f0 = 1e-3;      % Set the domain from [f0, fmax]
  8. fmax = 50;
  9.  
  10. % Leading-order solution 
  11. q0 = @(f) (f./(f+1)).^(1/2);
  12.  
  13. % Define the differential equation 
  14. dq = @(f, q) (-1i/2)*(q.^2 - q0(f).^2)./(ep*q^2.*q0(f));
  15.  
  16.  
  17. [phi, Q]  = ode45(dq, [f0, fmax], q0(f0));
  18. qsol = Q(:,1);
  19.  
  20. figure(1); 
  21. plot(phi, real(qsol), 'b');
  22. hold on
  23. plot(phi, q0(phi), 'r');
  24. hold off
  25. xlabel('\phi'); ylabel('q(\phi)');
  26. title('Solution of the model stern-wave problem');