I have been working on recreating the numerical solutions to the Moore-Saffman trailing vortex found in “Axial flow in laminar trailing vortices” by Moore and Saffman (1973) and “Linear stability of the Moore-Saffman model for a trailing wingtip vortex” and have been having an issue generating accurate results for the axial velocity profile (W_n) that they plot in the articles.

The issue is presenting itself after I use ode45 to solve for the particular solution to the following second order ODE.
enter image description here

Please take a look at my code let me know if you can help.

% Constant Declaration
xc = 0;
yc = 0;
B = .5; 
ufree = 1;
vis = 0.25;
z = 1;
time = z/ufree; % Assumption is made for small angle of attacks: Appendix A (Moore Saffman 1973)

delta = 0.005; % Distance between grid points
span = 6;
size = -span:delta:span;
numnodes = (span)/delta; % numbers from -5 to 5 time delta

r = 0.0000001:delta:span;
eta = zeros(1, length(r)); 
Vn = zeros(5, length(r));
Pn0 = zeros(5, length(r));
PnPrime = zeros(5, length(r));
Wn = zeros(5, length(r));

etaFun = @(y) -(y.^2)/(4*vis*time);

eta = arrayfun(etaFun, r);
% Solving Vn, Pn, and Wn for 5 different n values

loopvar = 1;

for j=1:loopvar
    % solving for 5 different n values
    % n = 0.2 + 0.1*j;
    n = 0.5;

    % solving Vn
    VnFun = @(x) (2.^(-n)) .* gamma((3/2)-(n/2)) .* ((-x).^(1/2)) .* hypergeom((1/2)+(n/2),2,x);
    Vn(j,:) = arrayfun(VnFun,eta);

    % solving Pn
    integrand = @(z, VnFun) (z.^(-1)) .* VnFun(z).^2;
    Pn0 = @(VnFun) -(1/2) * integral(@(z) integrand(z, VnFun), -Inf, 0);
    PnFun = @(x, VnFun) Pn0(VnFun) - (1/2) * integral(@(z) integrand(z, VnFun), 0, x);

    WnFun = @(eta, y, PnFun, integrand, VnFun) [y(2); -((1 - eta) * y(2) - n * y(1) + n * PnFun(eta, VnFun) + eta * integrand(eta, VnFun)) / eta];
    odeFun = @(eta, y) WnFun(eta, y, PnFun, integrand, @(x) VnFun(x));
    Wn0 = [0; -n*Pn0(VnFun)];
    
    etaCol = eta.';
    [t, WnTemp] = ode45(odeFun, etaCol, Wn0);

    WnI = WnTemp(:,1).';
    WnIslope = WnTemp(:,2).';
    WnAsymVar = ((2^(-1-2*n)) * ((1/n)-1));
    WnAsym = WnAsymVar .* (-eta).^(-n);

    intWnI = 0;
    % Area under WnI for finding Omega
    for m = numnodes-500:numnodes
        intWnI = intWnI + ((WnI(1,m-1)+WnI(1,m))/2) * (eta(1,m-1)-eta(1,m)) ;
    end
    
    sumWnIprime = 0;
    for m = numnodes-300:numnodes
        sumWnIprime = sumWnIprime + WnIslope(1,m);
    end
    dWnI = sumWnIprime/300;

    intOmega = ((-eta(1,numnodes-500))^(1-n) / (n-1)) - ((-eta(1,numnodes))^(1-n) / (n-1));
    WnIAsym = @(omega,eta) omega .* (-eta).^(-n); 

    %WnI and eta Trunction
    for m = 1:100
        WnItruncate(m) = WnI(numnodes-100+m);
        etatruncate(m) = eta(numnodes-100+m);
    end

    fit_data = [WnItruncate];


    options = optimoptions('lsqcurvefit', 'Display', 'iter');

    % omega least squares fit 
    [fit_params, res] = lsqcurvefit(WnIAsym, 1, etatruncate, fit_data, [], [], options);
    omega = fit_params;
    WnIAsymdata = WnIAsym(omega,eta);

    epsilonN = (WnAsymVar - omega)/gamma(1-n);
    
    for i = 1:numnodes
        Wn(j,i) = WnI(j,i) +  epsilonN * hypergeom(n,1,eta(1,i));
    end
    Wbar(j) = max(Wn(j,:)) - min(Wn(j,:));

    % vortex strength parameter
    % q(j) = ((vis * t)^(n/2) * ufree)/(2 * B * Wbar(j));
    % V(j,:) = 2 .* q .* Vn(j,:);
    % W(j,:) = (1/Wbar) .* Wn(j,:);

    % u_tangent(j,:) = (B ./ (vis * t).^(n/2) ) .* Vn(j,:);
    % u_axial(j,:) = B^2 ./ (ufree * (vis * t).^n) .* Wn(j,:);
end

I thought it was an issue with passing functions to ode45, but I made adjustments and still haven’t been able to find the solution. I can calculate correct azimuthal velocity profiles and it seems that Pn is producing correct values as well. I am expecting results that look like the following.
Numerical solution from Moore and Saffman (2014)

Here are the results I’ve been generating:
Current Azimuthal and Axial Velcity Profiles