function [T,Y] = TDRK4(f1,f2,td,a,b,ya,hist,n)
% DDE solver using 4th order Runge Kutta

h = (b - a)/n;

hstep = td/h;

T = linspace(a,b,n+1);

dim = size(ya,1);

hT = linspace(a-td,a-h,hstep);

hY = hist;

Y = zeros(dim,n+1);

Y(:,1) = ya;

for j=1:n
    if j-hstep < 1
        k1=feval(f1,T(j),Y(:,j)) + feval(f2,hT(j),hY(:,j));
        k2=feval(f1,T(j)+(h/2),Y(:,j)+(h/2)*k1) + feval(f2,hT(j)+(h/2),hY(:,j)+(h/2)*k1);
        k3=feval(f1,T(j)+(h/2),Y(:,j)+(h/2)*k2) + feval(f2,hT(j)+(h/2),hY(:,j)+(h/2)*k2);
        k4=feval(f1,T(j+1),Y(:,j)+h*k3) + feval(f2,hT(j)+h,hY(:,j)+h*k3);     
        
    else
        k1=feval(f1,T(j),Y(:,j)) + feval(f2,T(j-hstep),Y(:,j-hstep));
        k2=feval(f1,T(j)+(h/2),Y(:,j)+(h/2)*k1) + feval(f2,T(j-hstep)+(h/2),Y(:,j-hstep)+(h/2)*k1);
        k3=feval(f1,T(j)+(h/2),Y(:,j)+(h/2)*k2) + feval(f2,T(j-hstep)+(h/2),Y(:,j-hstep)+(h/2)*k2);
        k4=feval(f1,T(j+1),Y(:,j)+h*k3) + feval(f2,T(j+1-hstep),Y(:,j-hstep)+h*k3);

    end
    
    Y(:,j+1) = Y(:,j) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    
end

end

