I want to solve ODE for a complex systems in a rk4 method

I know to solve for simple system But i want to solve complex systems using matrix , array in rk4 method. I have a simple rk4 code , now i want to extend to a complex system. You can assume variables are vector .

function [x,y] =  rk4th(dydx,xo,xf,yo,h)

x = xo:h:xf ;
y = zeros(1,length(x));
y(1)= yo ;

for i = 1:(length(x)-1)
    k_1 = dydx(x(i),y(i));
    k_2 = dydx(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = dydx((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = dydx((x(i)+h),(y(i)+k_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
dydx = @(x,y) 3.*exp(-x)-0.4*y;
%[x,y] =  rk4th(dydx,0,100,-0.5,0.5);
%plot(x,y,'o-');

end

Answer 1

[James Tursa]

Just turn all the y( ) stuff into row or column vectors. E.g., suppose we use column vectors, then something like this ...

y = zeros(numel(yo),length(x));
y(:,1) = yo;
for i = 1:(length(x)-1)
    k_1 = dydx( x(i)      , y(:,i)           );
    k_2 = dydx( x(i)+0.5*h, y(:,i)+0.5*h*k_1 );
    k_3 = dydx( x(i)+0.5*h, y(:,i)+0.5*h*k_2 );
    k_4 = dydx( x(i)+    h, y(:,i)+    h*k_3 );
    y(:,i+1) = y(:,i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end

[mks]

how can I write the two different system of equation below here;

mu=0.0001;h=rand() * (0.3 - 0.15) + 0.15;

p=0.5;

q=[];q1=[];

c1=[];c2=[];

for pp=1:1

filename=append('AA',int2str(pp),'.mat');

load(filename)

B=A1;

[n m]=size(B);

for i=1:n

for j=1:m

if B(i,j)>0

B(i,j)=1;

else B(i,j)=0;

end

end

end

b=eye(m);

b1=eye(n);

b=eye(m);

b1=eye(n);

for i=1:m

for j=1:m

if i==j

b(i,j)=1;

else

b(i,j)=0.0;

end

end

end

for i=1:n

for j=1:n

if i==j

b1(i,j)=1;

else

b1(i,j)=0.0;

end

end

end

k1=sum(B,1);

k2=sum(B,2);

g=rand() * (1.2 - 0.8) + 0.8;

for ii=1:m

B1(:,ii)=(B(:,ii)./(k1(ii)^p))*g;

end

for ii=1:n

B2(ii,:)=(B(ii,:)./(k2(ii)^p))*g;

end

A= diag(rand(1,n)*(1.1-0.8)+0.8) + (rand(n)-eye(n))*(0.05-0.01) + 0.01*eye(n);

% ts=0:0.05:3;

y0=[];

y0=[rand(m,1); rand(n,1)];

y0=y0';

y0=reshape(y0,[1,m+n]);

p0=y0(1:m);

q0=y0(m+1:m+n);

x=p0; % initial pollinator abundance

y=q0; % initial plant abundance

z=0.0001*rand(m,1)'; %initial norm abundance

k3=0.18;

d=0.5;

m3=0.5;

dA=0.6;

muA=0.0001;

muP=0.0001;

l=0.14;

t=0;

t_max =500; dt=0.01;

m1=t_max/dt;

k=0;

k_max=1;

k_n=10;

dk=(k_max-k)/k_n;

%

for jj=1:k_n+1

t=0.0;

while t<t_max

for j=1:m

% for i=1:n

c1(j)=B1(:,j)'*y';

c1(j)=c1(j)/(1+h*c1(j)); % growth due to mutualism for pollinators

end

for j=1:n

c2(j)=B2(j,:)*x';

c2(j)=c2(j)/(1+h*c2(j)); % growth due to mutualism for plants

end

a1= rand() * (0.35 - 0.05) + 0.05;

a2= rand() * (0.35 - 0.05) + 0.05;

B3=A*(x'.*x); % competition faced by pollinators

B4=A*(y'.*y); % competition faced by plants

x=x+ (a1.*x-diag(B3)'+muA+c1.*x)*dt;

y=y+(a2.*y-k.*y-diag(B4)'+muP+c2.*y)*dt;

t=t+dt;

end

q=[q; k x y ];

k=k+dk;

save(append('data_norm_opt_pol',int2str(pp)),'q','-v7.3');

end

end

figure

% plot(q(:,1),q(:,2:26));

plot(q(:,1),q(:,27:51));

% hold on

% plot(q(:,1),q(:,52:76));

PLease complete the RK4 method to integration;