function [tTrace,y1,y2,yana] = nldamp(t0,y0,dt,nSteps)

% use: [tTrace,y1,y2,yana] = nldamp(0,[0 1],0.01,6280);
%      plot(tTrace,y1-yana)

y1 = zeros(1,nSteps);
y2 = zeros(1,nSteps);
tTrace = zeros(1,nSteps);

t = t0;
y = y0;

for k=1:nSteps
    
    k1 = rkFunc(t,y);
    k2 = rkFunc(t+0.5*dt,y+0.5*dt*k1);
    k3 = rkFunc(t+0.5*dt,y+0.5*dt*k2);
    k4 = rkFunc(t+dt,y+dt*k3);
    
    y = y + dt/6*(k1+2*k2+2*k3+k4);
    
    t = t+dt;

	tTrace(k) = t;
	y1(k) = y(1);
    y2(k) = y(2);
end
   
% add initial condition to vector

y1 = [y0(1) y1];
y2 = [y0(2) y2];
tTrace = [t0 tTrace];

ep = 0.1;

yana = 1./sqrt(0.75*ep*tTrace+1).*sin(tTrace);

end

function yPrime = rkFunc(t,y)

ep = 0.1;

yPrime = zeros(1,length(y));

yPrime(1) =  y(2);
yPrime(2) = -y(1)-ep*y(2)^3;
end