Why is my power series method not giving the correct solution?

I'm trying to solve the equation of motion for an object falling freely in a Laplacian atmosphere. My first method is an RK2, but I also wanted to solve the equation in a different way, using the exact equation as discussed in this paper: https://ttu-ir.tdl.org/server/api/core/bitstreams/e7af3650-0211-4a8e-aba8-862a7b755662/content.

The RK2 method seems to work just fine, but for some reason the other one doesn't. I think that there might be a problem with the way I implemented the infinite sum in my code (line 38-40, for k=1:50...).

Thank you for your help in advance!

This is my code:

function FreeFall(d, z0, h, m)
format long
T = 300; % interval length
N = (T-1)/h; % stepsize 
t = linspace(0, T, N+1); % grid
v = zeros(1, N+1); % vector of the numerical approximation
v_series = zeros(1, N+1); % vector of the power series approximation
dz = zeros(1, N+1); % vector for the vertical coordinate of the falling object
z = zeros(1, N+1); % vector for the vertical coordinate of the falling object (distance from the ground)
dz_series = zeros(1, N+1); % vector for the vertical coordinate of the falling object
z_series = zeros(1, N+1); % vector for the vertical coordinate of the falling object (distance from the ground)
q = zeros(1, 50); % vector for the sum of the power series
lambda = 7462.1; % characteristic height of the atmosphere [m]
g = 9.81; % gravity acceleration [m/s^2]
k0 = 0.22*(d^2); % coefficient of air resistance at sea level [m^2]
vt = sqrt((m*g)/k0); % terminal speed [m/s]
a = (2*g*lambda)/vt^2; % dimensionless constant
%%% computing the numerical approximation - RK2 method %%%
v(1) = 0; % initial condition
for j = 1:N
    dz(j) = v(j)*h; % calculating the height of the falling object
    z(j) = z0 - sum(dz(1:j));
    % RK2
    k1 = h*f(v(j), z(j), g, k0, lambda, m);
    k2 = h*f(v(j) + (k1/2), z0 - sum(dz(1:j)+((h^2*k1)/2)), g, k0, lambda, m);
    v(j+1) = v(j) + k2;
end
%%% approximation using the power series method %%%
v_series(1) = 0;
for i = 1:N
    dz_series(i) = v_series(i)*h;
    z_series(i) = z0 - sum(dz_series(1:i));
    for k = 1:50
        q(k) = (a^k * (exp(-(k*z_series(i))/lambda - exp(-(k*z0)/lambda)))) / (k*factorial(k));
    end
    v_series(i+1) = ( sqrt(a) * exp(-(a/2) * exp(-z_series(i)/lambda)) * sqrt((z0-z_series(i))/lambda + sum(q)) )*vt;
end

%%% plotting %%%
hfig = figure;
%plot(t, v, 'LineWidth', 1.5)
plot(t, v_series)
hold on;
plot(t, v)
%plot(1-(z/z0), v/vt)
%xlim([0, 1]); % set the plotting limit on the x-axis
%xlabel('1-z/z_{0}');
%ylabel('v/v_{t}');
xlabel('t [s]');
ylabel('v [m/s]');
%title('Speed as a function of altitude');
% fname='myfigure1';
% 
% picturewidth = 20; % set this parameter and keep it forever
% hw_ratio = 0.65; % feel free to play with this ratio
% set(findall(hfig,'-property','FontSize'),'FontSize',17) % adjust fontsize to your document
% 
% set(findall(hfig,'-property','Box'),'Box','off') % optional
% set(findall(hfig,'-property','Interpreter'),'Interpreter','latex') 
% set(findall(hfig,'-property','TickLabelInterpreter'),'TickLabelInterpreter','latex')
% set(hfig,'Units','centimeters','Position',[3 3 picturewidth hw_ratio*picturewidth])
% pos = get(hfig,'Position');
% set(hfig,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
% print(hfig,fname,'-dpdf','-painters','-fillpage')
%print(hfig,fname,'-dpng','-painters')
%%% 'f' function %%%
    function output = f(v, z, g, k0, lambda, m)
        output = -(-g + ( (k0*(v^2)*exp(-z/lambda)) / m));
    end
end

[William Rose]

@Vilmos Laszlo Balazs,

Please explain what fails when you call this function.

Please provide an example script, or a line of code, that calls the function.

When one defines a function, it is a helpful to include comments in the function that define each of the input parameters: what each parameter means, and the units for each.

What is the infinite power series you are trying to approximate with a finite sum? Does the error change when you change the number of terms in the finite sum approximation?

[Yash]

Hi Vilmos,

It seems you're facing issues while solving the mentioned equation of motion. To help resolve this, please consider the following points:

1. Provide a set of values for parameters "d", "z0", "h", and "m". Also provide their short description.
2. If the variable "N" represents the step size, which is not necessarily an integer (it's 299/h, depending on the value of "h"), it cannot be the size of the arrays declared below. Kindly make sure that you are using the correct "N".
3. Kindly provide a line of code that calls the function.

By addressing these points, we can work towards a better resolution.