help solve odefcn with time dependence

Hi,

Can you please help me solve my problem with this 16 coupled diffrential equation ? With Jkl-all time dependent series of size (1,991).

Thank you very much!

function dydt = odecdn(t,x,A,U,J,N)
    q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
    q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
    p1 = (abs(A(:,1))).^2;
    p2 = (abs(A(:,2))).^2;

    J11= transpose(real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));

    J12= transpose(real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));

    J13= transpose(real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));

    J14= transpose(real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));

    J21= transpose(real((-J/2).*sin(q2-q1)));

    J22= transpose(real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));

    J23= transpose(real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));

    J24= transpose(real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));

    J31= transpose(real(-J*sqrt(p1.*p2).*cos(q1-q2)));

    J32= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));

    J33= transpose(real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));

    J34= transpose(real(-0.5.*J*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));

    J41= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));

    J42= transpose(real(J*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));

    J43= transpose(real(J*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));

    J44= transpose(real(J*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
    dydt=zeros(9991,9991)
    y=zeros(16,9991)
    dydt(1)= (J11.*(y(1))+ (J12).*(y(5))+(J13).*(y(9))+(J14).*(y(13))); 
    dydt(2)= J11.*(y(2))+ (J12).*(y(6))+(J13).*(y(10))+(J14).*(y(14));
    dydt(3)= (J11.*(y(3))+ J12.*(y(7))+J13.*(y(11))+J14.*(y(15)));
    dydt(4)= (J11.*(y(4))+ J12.*(y(8))+J13.*(y(12))+J14.*(y(16)));
    dydt(5)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
    dydt(6)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
    dydt(7)= (J21.*(y(3))+ J22.*(y(7))+J23.*(y(11))+J24.*(y(15)));
    dydt(8)= (J21.*(y(4))+ J22.*(y(8))+J23.*(y(12))+J24.*(y(16)));
    dydt(9)= (J31.*(y(1))+ J32.*(y(5))+J33.*(y(9))+J34.*(y(13)));
    dydt(10)= (J31.*(y(2))+ J32.*(y(6))+J33.*(y(10))+J34.*(y(14)));
    dydt(11)= (J31.*(y(3))+ J32.*(y(7))+J33.*(y(11))+J34.*(y(15)));
    dydt(12)= (J31.*(y(4))+ J32.*(y(8))+J33.*(y(12))+J34.*(y(16)));
    dydt(13)= (J41.*(y(1))+ J42.*(y(5))+J43.*(y(9))+J44.*(y(13)));
    dydt(14)= (J41.*(y(2))+ J42.*(y(6))+J43.*(y(10))+J44.*(y(14)));
    dydt(15)= (J41.*(y(3))+ J42.*(y(7))+J43.*(y(11))+J44.*(y(15)));
    dydt(16)= (J41.*(y(4))+ J42.*(y(8))+J43.*(y(12))+J44.*(y(16)));
end

[Torsten]

The number of differential equations and the number of unknown functions must be the same.

Since you initialize

dydt=zeros(9991,9991)
y=zeros(16,9991)

you seem to have 16*9991 unknown functions and 9991*9991 differential equations.

This is not solvable.

Answer 1

[Torsten]

tspan = 1:0.1:10^3;
y0 = [sqrt(0.9) sqrt(1.1) sqrt(1)];
U=1;
J=1;
N=3;  % gewenste aantal deeltjes per site
mu=U*N-J; % chemische potentiaal
%lambda=J/UN
%k=J
%abs(y(length(t),1))^2+abs(y(length(t),2))^2
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
% grondtoestand vinden
[t,y] = ode15s(@(t,y) odefcnim(y,U,J,mu,N), tspan, y0, options);
%plot(t,abs(y(:,1:3)).^2)
% begintoestand
y0=y(end,:).*exp(1i*[1 0.5 -1].*0.1.*[sqrt(0.9) sqrt(1.1) sqrt(1)]);
% dynamica
[t,y] = ode15s(@(t,y) odefcn(y,U,J,N), tspan, y0, options);

%solve for A such that, dA/dt= JA , with J the jacobian matrix.
A=y;
%y=A=9991x3-matrix
tspan = [1:0.1:10^3];
ic = [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,s] = ode15s(@(t,y) odefcn2(t,y,tspan,A,U,J,N), tspan, ic, opts);
%odesolver calculates A in function of the time.


function dydt = odefcn(y,U,J,N)
dydt      = zeros(3,1);
dydt(1)   = -1i * (-J*0.5*y(2)+ U*N*abs(y(1))^2*y(1));
dydt(2)   = -1i * (-J*0.5*y(1)-J*0.5*y(3)+ U*N*abs(y(2))^2*y(2));
dydt(3)   = -1i * (-J*0.5*y(2)+U*N*(N-norm(y(1))^2-norm(y(2))^2)*y(3));
end
function dydt = odefcnim(y,U,J,mu,N)
dydt      = zeros(3,1);
dydt(1)   = -(-J*0.5*y(2)+ U*N*(abs(y(1))^2)*y(1) - mu*y(1)   );
dydt(2)   = -(-J*0.5*y(1)-J*0.5*y(3)+ U*N*(abs(y(2))^2)*y(2) - mu*y(2) );
dydt(3)   = -(-J*0.5*y(2)+U*N*((N-norm(y(1))^2-norm(y(2))^2))*y(3) - mu*y(3) );
end
function dydt = odefcn2(t,y,tspan,Amat,U,J,N)
A = interp1(tspan,Amat,t);
q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
p1 = (abs(A(:,1))).^2;
p2 = (abs(A(:,2))).^2;
J11= transpose(real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
J12= transpose(real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
J13= transpose(real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
J14= transpose(real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
J21= transpose(real((-J/2).*sin(q2-q1)));
J22= transpose(real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
J23= transpose(real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
J24= transpose(real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
J31= transpose(real(-J*sqrt(p1.*p2).*cos(q1-q2)));
J32= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));
J33= transpose(real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
J34= transpose(real(-0.5.*J*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
J41= transpose(real(J.*sqrt(p1.*p2).*cos(q1-q2)));
J42= transpose(real(J.*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
J43= transpose(real(J.*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
J44= transpose(real(J.*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
dydt = zeros(16,1);
dydt(1)= (J11.*(y(1))+ (J12).*(y(5))+(J13).*(y(9))+(J14).*(y(13))); 
dydt(2)= J11.*(y(2))+ (J12).*(y(6))+(J13).*(y(10))+(J14).*(y(14));
dydt(3)= (J11.*(y(3))+ J12.*(y(7))+J13.*(y(11))+J14.*(y(15)));
dydt(4)= (J11.*(y(4))+ J12.*(y(8))+J13.*(y(12))+J14.*(y(16)));
dydt(5)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(6)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(7)= (J21.*(y(3))+ J22.*(y(7))+J23.*(y(11))+J24.*(y(15)));
dydt(8)= (J21.*(y(4))+ J22.*(y(8))+J23.*(y(12))+J24.*(y(16)));
dydt(9)= (J31.*(y(1))+ J32.*(y(5))+J33.*(y(9))+J34.*(y(13)));
dydt(10)= (J31.*(y(2))+ J32.*(y(6))+J33.*(y(10))+J34.*(y(14)));
dydt(11)= (J31.*(y(3))+ J32.*(y(7))+J33.*(y(11))+J34.*(y(15)));
dydt(12)= (J31.*(y(4))+ J32.*(y(8))+J33.*(y(12))+J34.*(y(16)));
dydt(13)= (J41.*(y(1))+ J42.*(y(5))+J43.*(y(9))+J44.*(y(13)));
dydt(14)= (J41.*(y(2))+ J42.*(y(6))+J43.*(y(10))+J44.*(y(14)));
dydt(15)= (J41.*(y(3))+ J42.*(y(7))+J43.*(y(11))+J44.*(y(15)));
dydt(16)= (J41.*(y(4))+ J42.*(y(8))+J43.*(y(12))+J44.*(y(16)));
end

Answer 2

[Hélène De Witte]

Hi, I have tried dydt=zeros(9991,9991) y=zeros(16,9991) , tol, but even then it blocked. But the Jkl's are time dependent is there an option to get a counter,I have to solve: dM/dt = J*M and J and M are 4x4 matrices, and J is time fependent,I only know matrix J in function of the time-step. So I wrote the 16 deffrential equations but when I run it it blocks. Can you please help me a second time, please?

I'm really thankful! Thank you very much!

%solve for A such that, dA/dt= JA , with J the jacobian matrix.
A=y;
%y=A=9991x3-matrix
tspan = [1:0.1:10^3];
ic = [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ];
opts = odeset('RelTol',1e-20,'AbsTol',1e-20);
[t,s] = ode45(@(t,y) odefcn(t,y,A,U,J,N), tspan, ic, opts);
%odesolver calculates A in function of the time.

function dydt = odefcn(t,x,A,U,J,N)
    q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
    q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
    p1 = (abs(A(:,1))).^2;
    p2 = (abs(A(:,2))).^2;
    J11= transpose(real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
    J12= transpose(real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
    J13= transpose(real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
    J14= transpose(real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
    J21= transpose(real((-J/2).*sin(q2-q1)));
    J22= transpose(real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
    J23= transpose(real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
    J24= transpose(real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
    J31= transpose(real(-J*sqrt(p1.*p2).*cos(q1-q2)));
    J32= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));
    J33= transpose(real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
    J34= transpose(real(-0.5.*J*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
    J41= transpose(real(J.*sqrt(p1.*p2).*cos(q1-q2)));
    J42= transpose(real(J.*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
    J43= transpose(real(J.*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
    J44= transpose(real(J.*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
    dydt=zeros(16,9991)
    y=zeros(16,9991)
    dydt(1)= (J11.*(y(1))+ (J12).*(y(5))+(J13).*(y(9))+(J14).*(y(13))); 
    dydt(2)= J11.*(y(2))+ (J12).*(y(6))+(J13).*(y(10))+(J14).*(y(14));
    dydt(3)= (J11.*(y(3))+ J12.*(y(7))+J13.*(y(11))+J14.*(y(15)));
    dydt(4)= (J11.*(y(4))+ J12.*(y(8))+J13.*(y(12))+J14.*(y(16)));
    dydt(5)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
    dydt(6)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
    dydt(7)= (J21.*(y(3))+ J22.*(y(7))+J23.*(y(11))+J24.*(y(15)));
    dydt(8)= (J21.*(y(4))+ J22.*(y(8))+J23.*(y(12))+J24.*(y(16)));
    dydt(9)= (J31.*(y(1))+ J32.*(y(5))+J33.*(y(9))+J34.*(y(13)));
    dydt(10)= (J31.*(y(2))+ J32.*(y(6))+J33.*(y(10))+J34.*(y(14)));
    dydt(11)= (J31.*(y(3))+ J32.*(y(7))+J33.*(y(11))+J34.*(y(15)));
    dydt(12)= (J31.*(y(4))+ J32.*(y(8))+J33.*(y(12))+J34.*(y(16)));
    dydt(13)= (J41.*(y(1))+ J42.*(y(5))+J43.*(y(9))+J44.*(y(13)));
    dydt(14)= (J41.*(y(2))+ J42.*(y(6))+J43.*(y(10))+J44.*(y(14)));
    dydt(15)= (J41.*(y(3))+ J42.*(y(7))+J43.*(y(11))+J44.*(y(15)));
    dydt(16)= (J41.*(y(4))+ J42.*(y(8))+J43.*(y(12))+J44.*(y(16)));
end

[Torsten]

What are the time instants that correspond to the 9991 entries of your matrix A ?

[Hélène De Witte]

Thank you very much! Well, it's also a solution (I try to solve a biger problem and needed the y first, it were also diffrential equations) of a ode function defined as follow:

tspan = 1:0.1:10^3;
y0 = [sqrt(0.9) sqrt(1.1) sqrt(1)];
U=1;
J=1;
N=3;  % gewenste aantal deeltjes per site

mu=U*N-J; % chemische potentiaal
%lambda=J/UN
%k=J
%abs(y(length(t),1))^2+abs(y(length(t),2))^2
options = odeset('RelTol',1e-10,'AbsTol',1e-10);

% grondtoestand vinden
[t,y] = ode15s(@(t,y) odefcnim(y,U,J,mu,N), tspan, y0, options);

%plot(t,abs(y(:,1:3)).^2)

% begintoestand
y0=y(end,:).*exp(1i*[1 0.5 -1].*0.1.*[sqrt(0.9) sqrt(1.1) sqrt(1)]);
% dynamica
[t,y] = ode15s(@(t,y) odefcn(y,U,J,N), tspan, y0, options);

function dydt = odefcn(y,U,J,N)
  dydt      = zeros(3,1);
  dydt(1)   = -1i * (-J*0.5*y(2)+ U*N*abs(y(1))^2*y(1));
  dydt(2)   = -1i * (-J*0.5*y(1)-J*0.5*y(3)+ U*N*abs(y(2))^2*y(2));
  dydt(3)   = -1i * (-J*0.5*y(2)+U*N*(N-norm(y(1))^2-norm(y(2))^2)*y(3));
end

function dydt = odefcnim(y,U,J,mu,N)
  dydt      = zeros(3,1);
  dydt(1)   = -(-J*0.5*y(2)+ U*N*(abs(y(1))^2)*y(1) - mu*y(1)   );
  dydt(2)   = -(-J*0.5*y(1)-J*0.5*y(3)+ U*N*(abs(y(2))^2)*y(2) - mu*y(2) );
  dydt(3)   = -(-J*0.5*y(2)+U*N*((N-norm(y(1))^2-norm(y(2))^2))*y(3) - mu*y(3) );
end

Answer 3

[Hélène De Witte]

Thank you very much! Maybe this helps too, In fact, I can solve it, like this, but I want to understand how I can fix it without manupulating the problem, because it ask a lot of time...

s=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]
for stap=1:2
    teller=1
    tspan = [teller:0.1:teller+0.1];
    ic = [s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end)];
    opts = odeset('RelTol',1e-20,'AbsTol',1e-20);
    [t,s] = ode45(@(t,y) kale(t,y,A,U,J,N,stap), tspan, ic, opts);
    %odesolver calculates A in function of the time.
    teller=teller+0.1
end

function dydt = kale(t,y,A,U,J,N, stap)
    q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
    q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
    p1 = (abs(A(:,1))).^2;
    p2 = (abs(A(:,2))).^2;
    J11= (real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
    J12= (real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
    J13= (real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
    J14= (real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
    J21= (real((-J/2).*sin(q2-q1)));
    J22= (real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
    J23= (real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
    J24= (real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
    J31= (real(-J*sqrt(p1.*p2).*cos(q1-q2)));
    J32= (real(J.*sqrt(p1.*p2).*cos(q1-q2)));
    J33= (real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
    J34= (real(-0.5.*J.*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
    J41= (real(J.*sqrt(p1.*p2).*cos(q1-q2)));
    J42= (real(J.*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
    J43= (real(J.*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
    J44= (real(J.*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
    dydt= zeros(16,1);
    dydt(1)= (J11(stap,1)*(y(1))+ (J12(stap,1))*(y(5))+(J13(stap,1))*(y(9))+(J14(stap,1))*(y(13))); 
    dydt(2)= J11(stap,1)*(y(2))+ (J12(stap,1))*(y(6))+(J13(stap,1))*(y(10))+(J14(stap,1))*(y(14));
    dydt(3)= (J11(stap,1)*(y(3))+ J12(stap,1)*(y(7))+J13(stap,1)*(y(11))+J14(stap,1)*(y(15)));
    dydt(4)= (J11(stap,1)*(y(4))+ J12(stap,1)*(y(8))+J13(stap,1)*(y(12))+J14(stap,1)*(y(16)));
    dydt(5)= (J21(stap,1)*(y(1))+ J22(stap,1)*(y(5))+J23(stap,1)*(y(9))+J24(stap,1)*(y(13)));
    dydt(6)= (J21(stap,1)*(y(1))+ J22(stap,1)*(y(5))+J23(stap,1)*(y(9))+J24(stap,1)*(y(13)));
    dydt(7)= (J21(stap,1)*(y(3))+ J22(stap,1)*(y(7))+J23(stap,1)*(y(11))+J24(stap,1)*(y(15)));
    dydt(8)= (J21(stap,1)*(y(4))+ J22(stap,1)*(y(8))+J23(stap,1)*(y(12))+J24(stap,1)*(y(16)));
    dydt(9)= (J31(stap,1)*(y(1))+ J32(stap,1)*(y(5))+J33(stap,1)*(y(9))+J34(stap,1)*(y(13)));
    dydt(10)= (J31(stap,1)*(y(2))+ J32(stap,1)*(y(6))+J33(stap,1)*(y(10))+J34(stap,1)*(y(14)));
    dydt(11)= (J31(stap,1)*(y(3))+ J32(stap,1)*(y(7))+J33(stap,1)*(y(11))+J34(stap,1)*(y(15)));
    dydt(12)= (J31(stap,1)*(y(4))+ J32(stap,1)*(y(8))+J33(stap,1)*(y(12))+J34(stap,1)*(y(16)));
    dydt(13)= (J41(stap,1)*(y(1))+ J42(stap,1)*(y(5))+J43(stap,1)*(y(9))+J44(stap,1)*(y(13)));
    dydt(14)= (J41(stap,1)*(y(2))+ J42(stap,1)*(y(6))+J43(stap,1)*(y(10))+J44(stap,1)*(y(14)));
    dydt(15)= (J41(stap,1)*(y(3))+ J42(stap,1)*(y(7))+J43(stap,1)*(y(11))+J44(stap,1)*(y(15)));
    dydt(16)= (J41(stap,1)*(y(4))+ J42(stap,1)*(y(8))+J43(stap,1)*(y(12))+J44(stap,1)*(y(16)));
end

[Hélène De Witte]

Yes, this works!Thank you very much!