solving coupled system of odes

hello everyone i am trying to solve two second order odes using runge kutta 4 and adams bashforth4 method. I'm not sure what part im doing wrong or how i should change it. The two odes in the system are: dv1/dt=-3x1+2x2-0.05v1+0.02 v2, dx1/dt=v1 and dv2/dt=0.8x1-2x2+0.008v1-0.01v2

Answer 1

[Torsten]

Try this (you can run it as is):

% this is the main function in another file :
function main
%tspan = [0 500];
%y0 = [-2;1;5;-3];
%[t,y] = ode15s(@f_MYODE,tspan,y0);
%figure(1)
%plot(t,y(:,1))
t0=0;   
tf=500;
z0(1)=-2;
z0(2)=1;
z0(3)=5;
z0(4)=-3;
NS=2000;
RHS=@f_MYODE; 
[yRK, tRK] = RK4 ( RHS, t0, z0, tf, NS);
plot(tRK,yRK(1,:))
NS=20000;
[yAB,tAB]  = AB4(RHS, t0, z0, tf, NS);
hold on
plot(tAB,yAB(1,:))
%figure1 = figure;
% Create axes
%axes1 = axes('Parent',figure1,'FontSize',20,'FontName','Times New Roman'); 
%box(axes1,'on');
%grid(axes1,'on');
%hold(axes1,'all');
% Create xlabel
%xlabel('x','FontSize',20,'FontName','Times New Roman'); 
% Create ylabel
%ylabel('y','FontSize',20,'FontName','Times New Roman');
% Create multiple lines using matrix input to plot
%plot1 = plot(tRK,yRK(1,:),'Parent',axes1,'LineWidth',3);
%plot2 = plot(tAB,yAB(1,:),'Parent',axes1,'LineWidth',3);
%set(plot1(1),'Marker','o','Color',[1 0 0],'DisplayName','RK4');
%set(plot2(1),'Marker','d','Color',[0 1 0],'DisplayName','AB4');
% Create legend
%legend(axes1,'show');
end
%and this is the dertivative of the function in another file 
function f=f_MYODE(t,z) %Derivative of the problem
f=[z(2);-3.*z(1)+2.*z(3)-0.05.*z(2)+0.02.*z(4);z(4);0.8.*z(1)-2.*z(3)+0.008.*z(2)-0.012.*z(4)];
end
% this is for rk4, should be saved in one file
function [wi, ti] = RK4 ( RHS, t0, x0, tf, N )
if size(x0,2)>1
    x0=x0';
end
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = zeros( neqn, N+1 );
wi(:, 1) = x0;
h = ( tf - t0 ) / N;
for i = 1:N
    k1 = RHS( ti(i) , wi(:,i));
    k2 = RHS( ti(i) + h/2 , wi(:,i) + k1*h/2);
    k3 = RHS( ti(i) + h/2 , wi(:,i) + k2*h/2);
    k4 = RHS( ti(i) + h , wi(:,i) + k3*h);
    wi(:,i+1) = wi(:,i) + h*( k1 + 2*k2 + 2*k3 + k4)/6;
end
end
%this is ab4 saved in another file
function [wi, ti] = AB4 ( RHS, t0, x0, tf, N )
if size(x0,2)>1
    x0=x0';
end
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = zeros( neqn, N+1 );
wi(:, 1) = x0;
h = ( tf - t0 ) / N;
save_AB = zeros(neqn,3);
for i = 1:3
    k1 = RHS( ti(i) , wi(:,i));
    k2 = RHS( ti(i) + h/2 , wi(:,i) + k1*h/2);
    k3 = RHS( ti(i) + h/2 , wi(:,i) + k2*h/2);
    k4 = RHS( ti(i) + h , wi(:,i) + k3*h);
    wi(:,i+1) = wi(:,i) + h*( k1 + 2*k2 + 2*k3 + k4)/6;
    save_AB(:,i) = k1;
end
for i=4:N
    save_AB_new = RHS(ti(i),wi(:,i));
    wi(:,i+1) = wi(:,i) + h*( 55*save_AB_new - 59*save_AB(:,3) + ...
                              37*save_AB(:,2) - 9*save_AB(:,1))/24;
    save_AB(:,1) = save_AB(:,2);
    save_AB(:,2) = save_AB(:,3);
    save_AB(:,3) = save_AB_new;
end    
end

[Jackie Longman]

Thank you for your answer, my main function works now, but I dont know how I should be getting x1 and x2, could you help me with that? What I mean is that I assume the answer that Im getting is x1 but how do I change my code to get x2

[Torsten]

So you didn't know that wi(i,:) is z(i) ? But you coded it this way ... I'm really surprised.

[Jackie Longman]

Well it is because my tutor gave us the code for the other question and asked us to change it for a second order ode without explaining, thats why Im so confused. So I should just be changing it to (3,:) is the plotting part?

[Sam Chak]

Sometimes the new staff will inherit the codes from the manager, the professor, a scientist, or colleagues, and then tasked to modify some parameters and the structure of the codes to suit the needs of the project. Most of the time, the codes do not come with comments, making tracing the variables difficult.

[Torsten]

@Jackie Longman

Yes, if you want to plot z(3), change the 1 to 3 in the plotting part.

[Jackie Longman]

Thank you so much ! much appriciated