Matlab - Por favor ayuda con metodo de interpolación dinámica

%% 4

clc, clear, close all

m=input('Ingrese masa dividida por la aceleración de la gravedad: ');

k=input('Ingrese rigidez: ');

wn=input('Ingrese frecuencia natural (rad/s): ');

dzeta=input('Ingrese razon del amortiguamiento menor al critico en decimal: ');

deltat=input('Ingrese tamaño del paso del tiempo (seg): ');

uo=0 %Desplazamiento inicial

uuo=0 %Velocidad inicial

t=0:deltat:1 %Caso partícular ejemplo 5.1 Chopra

n=length(t)

pt=10*sin(pi*t/0.6) %Caso partícular ejemplo 5.1 Chopra

pt(pt<0)=0

%Cálculos iniciales

q=exp(-dzeta*wn*deltat)

wD=wn*sqrt(1-dzeta^2)

r=sin(wD*deltat)

p=cos(wD*deltat)

%COEFICIENTES EN LAS FÓRMULAS DE RECURRENCIA (? < 1)

A=q*((dzeta/sqrt(1-dzeta^2))*r+p)

B=q*((1/wD)*r)

C=(1/k)*(((2*dzeta)/(wn*deltat))+q*((((1-2*dzeta^2)/(wD*deltat))-(dzeta/sqrt(1-dzeta^2)))*r-(1+(2*dzeta)/(wn*deltat))*p))

D=(1/k)*(1-((2*dzeta)/(wn*deltat))+q*((((2*dzeta^2-1)/(wD*deltat)))*r+((2*dzeta)/(wn*deltat))*p))

Ap=-q*((wn/(sqrt(1-dzeta^2)))*r)

Bp=q*(p-(dzeta/(sqrt(1-dzeta^2)))*r)

Cp=(1/k)*((-1/deltat)+q*(((wn/(sqrt(1-dzeta^2)))+dzeta/(deltat*sqrt(1-dzeta^2)))*r+(1/deltat)*p))

Dp=(1/(k*deltat))*(1-q*((dzeta/(sqrt(1-dzeta^2)))*r+p))

Cpi=pt*C

Dpi=pt(2:n)*D

Cppi=pt*Cp

Dppi=pt(2:n)*Dp

i=2

for k=1

    Aui=A*uo

    Apui=Ap*uo

    Buui=B*uuo

    Bpuui=Bp*uuo

    for k=2:n

        ui=Aui+Buui+Cpi(k-1)+Dpi(k-1)

        uui=Apui+Bpuui+Cppi(k-1)+Dppi(k-1)

        Aui=A*ui

        Apui=Ap*ui

        Buui=B*uui

        Bpuui=Bp*uui

    end

end