How to code integral part in diff. eqn with 0 to inf limit ?

I want to code the above differential equation in matlab. But the only problem I am facing is how to code the integral part. I tried a lot but didn't get success. Please see below my trial code. Many thanks in advance.

function unbounded
clear;

global theta
t0 = 0; t1 = 5;
    options=ddeset('MaxStep',0.001,'RelTol',6,'AbsTol',7);
    sol = ddesd(@qvnn, @delay, [-0.2 0.5 0.4 -0.4 0 0 0 0 0.5 0.5 0.5 0.5 0 0 0 0 0.1 0.9], [t0, t1], options);
 %  plot(sol.y(17,:),sol.y(18,:),'r');
   plot(sol.y(17,:),'r');
    %xlabel('$q_1(t)$','Interpreter','latex');
 %ylabel('$q_2(t)$','Interpreter','latex');
%  legend('q_1(t)','q_2(t)','Location','northeast')
  % plot(sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'c',sol.x,sol.y(7,:),'r',sol.x,sol.y(8,:),'b');
function d = delay(t,~)
     d = [t-exp(t)./(exp(t)+1),5];
end
function dy = qvnn(t,y,Z)
         ylag = Z(:,1);
         
         ylag2 = Z(:,2);
         dy = zeros(18,1);
%first system
dy(1)=y(3);
dy(3)=-8*y(3)-(1+1/(y(1).*y(1)+1))*(0.5*y(1)-0.125*tanh(y(1))-0.175*tanh(y(2))-0.1*tanh(ylag(1))-0.15*tanh(ylag(2)))+0.2*y(5)+0.4*y(6);
dy(2)=y(4);
dy(4)=-10.4*y(4)-(1+1/(y(2).*y(2)+1))*(0.5*y(2)+0.4*tanh(y(1))+0.2*tanh(y(1))+0.05*tanh(y(2))+0.0755*tanh(ylag(1))+0.05*tanh(ylag(2)))-0.3*y(7)+0.6*y(8);


% second system
dy(9)=y(10);
dy(10)=-8*y(10)-(1+1/(y(9).*y(9)+1))*(0.5*y(9)-0.125*tanh(y(9))-0.175*tanh(y(11))-0.1*tanh(ylag(9))-0.15*tanh(ylag(11)))+0.2*y(13)+0.4*y(14)-10*y(17);
dy(11)=y(12);
dy(12)=-10.4*y(12)-(1+1/(y(11).*y(11)+1))*(0.5*y(11)+0.4*tanh(y(9))+0.2*tanh(y(9))+0.05*tanh(y(11))+0.0755*tanh(ylag(9))+0.05*tanh(ylag(11)))-0.3*y(15)+0.6*y(16)-12*y(18);


%%%%%%%%%%%%%%%% integral part %%%%%%%
a = y(1);
b = y(2);
c = y(9);
d = y(11);
f1 = @(t, theta)sin(2*theta)./(1+theta^2).*tanh(a(t-theta));
f2 = @(t, theta)sin(2*theta)./(1+theta^2).*tanh(b(t-theta));
f3 = @(t, theta)sin(2*theta)./(1+theta^2).*tanh(c(t-theta));
f4 = @(t, theta)sin(2*theta)./(1+theta^2).*tanh(d(t-theta));

 y(5) = integral(f1,0, inf);
 y(6) = integral(f2,0, inf);
 y(7) = integral(f1,0, inf);
 y(8) = integral(f2,0, inf);
 y(13) = integral(f3,0, inf);
 y(14) = integral(f4,0, inf);
 y(15) = integral(f3,0, inf);
 y(16) = integral(f4,0, inf);
 
 
 
%%%%%%%%%% we don't know how to define the unbounded intgral part which is give below %%%%%%%%%%%%%%
% y(5)=integration_{0,Inf}*sin(2*theta)./(1+theta^2)*tanh(tdiff(1))*dtheta;
% y(6)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(2)(t-\theta))* d\theta;
% y(7)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(1)(t-\theta))* d\theta;
% y(8)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(2)(t-\theta))* d\theta;
% y(13)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(9)(t-\theta))* d\theta;
% y(14)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(11)(t-\theta))* d\theta;
% y(15)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(9)(t-\theta))* d\theta;
% y(16)= integration_{0,Inf}*sin(2*theta)/(1+theta^2)*tanh(y(11)(t-\theta))* d\theta;
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%    error= y(t)-x(t-\tau)   %%%%%%%%%%%%%%%%%
dy(17)=dy(9)-dylag2(1);
dy(18)=dy(11)-dylag2(2);
end
end

[Torsten]

You must have access to the complete solution of the w_j(theta) values for theta <= t.

This is hard to achieve for MATLAB's ode solvers.

And why do you have 18 ODEs ? According to what you wrote above, you only have 4 since n=2.

[Rakesh]

Many thanks for your comments.

The integral part is a convolution type integration in which w() is an unknown function. w() is itself the solution of the differential equation so we don't know w_j(theta) values for theta <=t.

Actually, I am solving two similar differential equations (master and slave systems for synchronization) in the code that is why I wrote 18 ODEs. If somebody help me to find the solution of only the above diff. Eq then it is sufficient for my pupose.

Or

If I get information that how will I code the integral part then it will be sufficient for me.

[Torsten]

w() is itself the solution of the differential equation so we don't know w_j(theta) values for theta <=t.

We know w_j(theta) for theta <=t since we solve for w_j, and t is the time instant that we have reached during the solution process. But as I said: Accessing the solution w_j for the times smaller than the time t reached during the solution process and working with it (e.g. to calculate the convolution integral) is hard to archieve.

Since you must have access to the complete solution for values of theta <= t, using a delay differential equations solver is not necessary. Once you find a way to access the solution for all time values in the past, you can also use a standard ODE integrator (e.g. ODE45).

[Walter Roberson]

Why is Kj written that way, implying that they are different for different j, but they are defined independent of j? Why isn't the integral calculated once and factored out like integral K(theta) * sum of cij??

[Torsten]

Why isn't the integral calculated once and factored out like integral K(theta) * sum of cij??

The integrand is continued in the next line of the problem description - it's a kind of convolution integral.

[Rakesh]

@Walter Roberson Why is Kj written that way, implying that they are different for different j, but they are defined independent of j?

"j" is just a suffix denotes K_1, K_2,... but actually all are functions of \theta. In our case, we have taken K_1 = K_2.

Thank you @Torsten to answer Walter Roberson for his second question.

[Torsten]

Maybe you want to make a search whether in the meantime, there are solvers for integro-differential equations in the File Exchange.

[Rakesh]

@Torsten I am searching for such a solver in the File Exchange. Hoping to get the solution.

[Torsten]

To be honest: I doubt you will find a suitable solver since - additionally to the integral - your equation contains delays for the w_i.

If it were only the integral part, IDSOLVER might have worked:

https://www.sciencedirect.com/science/article/abs/pii/S0010465513003093?via%3Dihub