% Function Subprogram srung3.m
%
% Solving the nonlinear parts of the coupled differential
% equations using a fourth-order Runge-Kutta algorithm
%

function y = srung3(rho1,rho2,sig1,sig2,a,b,h,idum)

   ka1 = h * funa3(rho1, sig1, a, b);
   kb1 = h * funb3(rho2, sig2, a, b);
   ka2 = h * funa3(rho1, sig1, a+ka1/2, b+kb1/2);
   kb2 = h * funb3(rho2, sig2, a+ka1/2, b+kb1/2);
   ka3 = h * funa3(rho1, sig1, a+ka2/2, b+kb2/2);
   kb3 = h * funb3(rho2, sig2, a+ka2/2, b+kb2/2);
   ka4 = h * funa3(rho1, sig1, a+ka3, b+kb3);
   kb4 = h * funb3(rho2, sig2, a+ka3, b+kb3);
   
   if idum == 0
      y = a + (ka1+2*ka2+2*ka3+ka4)/6;
   else
      y = b + (kb1+2*kb2+2*kb3+kb4)/6;
end

%==============================================================
% Function Subprogram funa3.m
% Calculates the right-hand-side of the first part of
% the coupled differential equations.

function y = funa3(rho, sig, a, b)
y = -i * rho * conj(a) .* b - i * sig * ((abs(a)).^2 .* a + 2*(abs(b)).^2 .* a);

%==============================================================
% Function Subprogram funb3.m
% Calculates the right-hand-side of the second part of
% the coupled differential equations.
function y = funb3(rho, sig, a, b);
y = -i * rho * (1/2) * a.^2 - i * sig * ((abs(b)).^2 .* b + 2*(abs(a)).^2 .* b);