Matlab function to calculate exact vertical position and velocity of a single bouncing ball (with given coefficient of restitution less than one) at any time or array of times

ball_bounce.m

%@EliasHasle
%There probably is room for more optimization in the form of vectorization.
%You see, Matlab is not my first language.

function [ values ] = ball_bounce( initials, k, tt)
    %Input:
    %initials: column vector of initial vertical position and velocity
    %k: coefficient of restitution (must be less than one)
    %tt: array of times
    %
    %Output:
    %values: array of column vectors of vertical position and velocity at times tt

    g = 9.81; %I will convert to negative acceleration as needed.
    
    %initial conditions:
    yStart = initials(1);
    vStart = initials(2);

    %theoretical last bouncing point before (or at) start
    %(negative solution to constant-acceleration equation)
    t0 = (vStart - sqrt(vStart^2 + 2*g*yStart))/g;
    
    %Theoretical post-bounce velocity at t0.
    v0 = vStart-g*t0;
    
    %vn = v0*k^(n-1);
    
    %Time calculations:
    %dtn is the time between bounce n and bounce n+1
    %dtn = 2*vn/g = 2*(v0*k^n)/g = (2*v0/g)*k^n
    dt0 = (2*v0/g);
    %then dtn = dt0*k^n
    
    %tn = t0 + sum(dtn, 0, n-1) = t0 + dt0*sum(k^n, 0, n-1)
    % = t0 + dt0*(1-k^(n-1))/(1-k)
    tnfun = @(n) t0 + dt0*(1-k^(n-1))/(1-k);
    
    %Note that if k<1 (which it should be), this will converge to:
    T = t0 + dt0/(1-k);
    %That means the velocity will be exactly 0 for t>=T
    %(after infinite bouncing before T).
    
    %To find the right tn, I will need the inverse function n(t):
    %t - t0 = dt0*(1-k^(n(t)-1))/(1-k)
    %...
    nfun = @(t) floor(log(1-(t-t0)*(1-k)/dt0)/log(k)+1);
    
    %Allocate space for output:
    y = zeros(1,length(tt(1,:)));
    v = zeros(1,length(tt(1,:)));
        
    for i = 1:length(tt(1,:))
        t = tt(i);
        if (t>=T)
            y(1,i) = 0;
            v(1,i) = 0;
        else
            n = nfun(t);
            tn = tnfun(n);
            dt = t-tn;
            vn = v0*k^(n-1);
            y(1,i) = vn*dt - 0.5*g*dt^2;
            v(1,i) = vn - g*dt;
        end
    end
    
    values = [y;v];
end