MATLAB: Interp1 returning NaN for ode45

Question

The interp1 function returns a NaN value for the time-dependent functions that are later introduced in the ode45 function. Using pchip or spline instead of interp1 returns a value other than NaN, but the solutions are not valid.

My code is the one that follows, thank you in advance.

function [T_V , W_V , x_V , y_V , h_V , Vfin] = viraje(x0 , y0 , h0 , W0 ,...
    V0 , Vf , gamma0 , gammaf , chi0 , chif , chip ,  CD0 , k , S , cE)
clc
g = 9.80665;   %m/s2
tfin = (chif - chi0)/chip;
gammap = (gammaf-gamma0)/tfin;
Vp     = (Vf-V0)/tfin;
t_interp = linspace(0 , tfin , 10);
V        = V0     + Vp*t_interp;
gamma    = gamma0 + gammap*t_interp;
mu       = atan(chip/g*V.*cos(gamma)./(V/g*gammap + cos(gamma)));
h        = Hviraje(t_interp , h0 , V0 , Vp , gamma0 , gammap);
rho      = densidad(h);
      function dWdt = dif(t,W,t_interp,V,gamma,mu,rho)
          if Vp == 0
              V = V0;
          else
              V = interp1(V,t_interp,t,'pchip');
          end
          if gammap == 0
              gamma = gamma0;
          else
              gamma = interp1(gamma,t_interp,t,'pchip');
          end
          if Vp == 0 && gammap == 0
              mu = mu(1);
          else
              mu = interp1(mu,t_interp,t,'pchip');
          end
          if gamma0 == 0 && gammap == 0
              rho = densidad(h0);
          else
              rho = interp1(rho,t_interp,t,'pchip');
          end
          dWdt  = -cE*(0.5*rho.*(V.^2)*S*CD0 + 2*k./(rho*S)*(chip*W.*cos(gamma)./...
              (g*sin(mu))).^2 + W.*sin(gamma) + Vp/g*W)
      end
tspan = [0 tfin];
ic    = W0;
[t,W] = ode45(@(t,W) dif(t,W,t_interp,V,gamma,mu,rho) , tspan , ic);
W_V = W;
h_V = Hviraje(t , h0 , V0 , Vp , gamma0 , gammap);
x_V = Xviraje(t , x0 , V0 , Vp , gamma0 , gammap , chi0 , chip);
y_V = Yviraje(t , y0 , V0 , Vp , gamma0 , gammap , chi0 , chip);
rho   = densidad(h_V)';
gamma = gamma0 + gammap*t;
V     = V0 + Vp*t;
mu    = atan(chip/g*V.*cos(gamma)./(V/g*gammap + cos(gamma)));
T_V   = 0.5*rho.*V.^2*S*CD0 + 2*k./(rho*S).*(chip*W.*cos(gamma)./...
    (g*sin(mu))).^2 + W.*sin(gamma) + Vp/g*W;
Vfin = Vf;
end
function H = Hviraje(T , h0 , V0 , Vp , gamma0 , gammap)
      syms h(t)
      eqn     = diff(h,t) == (V0 + Vp*t)*sin(gamma0 + gammap*t);
      cond    = h(0) == h0;
      ySol(t) = dsolve(eqn,cond);
      H       = ySol(T);
end
function X = Xviraje(T , x0 , V0 , Vp , gamma0 , gammap , chi0 , chip)
      syms x(t)
      eqn     = diff(x,t) == (V0 + Vp*t)*cos(gamma0 + gammap*t)*cos(chi0 + chip*t);
      cond    = x(0) == x0;
      ySol(t) = dsolve(eqn,cond);
      X       = ySol(T);
end
function Y = Yviraje(T , y0 , V0 , Vp , gamma0 , gammap , chi0 , chip)
      syms y(t)
      eqn     = diff(y,t) == (V0 + Vp*t)*cos(gamma0 + gammap*t)*sin(chi0 + chip*t);
      cond    = y(0) == y0;
      ySol(t) = dsolve(eqn,cond);
      Y       = ySol(T);
end

Answer 1

V = interp1(V,t_interp,t,'pchip');

It looks arguments are swapped; I think it should be

V = interp1(t_interp,V,t,'pchip');

and similar for other INTERP1 calls.