2D ODE with constant? how to solve

Answer 1

[Sulaymon Eshkabilov]

Most parts of your code is ok, but within the loop, you have overlooked sth and thus, you final solutions are not quite accurate. Here is ODE45 simulation which can be compared with your simulation results.

ICs=[0.6;0.6];
a=0.10;
b=10;
t=[0,60];
F = @(t, z)([a-z(1)+z(1).^2*z(2);b-z(1).^2*z(2)]);
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-9);
[time, z]=ode45(F, t, ICs, OPTs);
figure(2)
plot(time,z(:,1),'b',time,z(:,2),'r')
xlabel('time')
ylabel('x(t) y(t)')
legend('x(t)', 'y(t)', 'location', 'best')
title('Schnackenberg eqn simulation'), xlim([0, 5])

figure(1)
plot(z(:,1),z(:,2),'k')
title('Simulation using ODE45'), grid on
xlabel('x(t)')
ylabel('y(t)')

Answer 2

[Sulaymon Eshkabilov]

Use odex (ode23, ode45, ode113, etc.) solvers. See this doc how to employ them in your exercise: https://www.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=srchtitle

[mays rashad]

Is this solution correct?

%x'=a-x+x^2y   y'=b-x^2y
clear all,close all, clc
x(1)=0.6;
y(1)=0.6;
a=0.10;
b=10;
h=0.02;
t=0:h:60;
for i=1:(length(t)-1)
   k1=h*(a-x(i)+y(i)*x(i)^2);
   L1=h*(b-y(i)*x(i)^2);
   k2=h*(a-(x(i)+k1/2)+(y(i)+L1/2)*(x(i)^2+k1/2));
   L2=h*(b-(y(i)+L1/2)*(x(i)^2+k1/2));
   k3=h*(a-(x(i)+k2/2)+(y(i)+L2/2)*(x(i)^2+k2/2));
   L3=h*(b-(y(i)+L2/2)*(x(i)^2+k2/2));
   k4=h*(a-(x(i)+k3)+(y(i)+L3)*(x(i)^2+k3));
   L4=h*(b-(y(i)+L3)*(x(i)^2+k3));
   x(i+1)=x(i)+(k1+2*k2+2*k3+k4)*(h/6);
   y(i+1)=y(i)+(L1+2*L2+2*L3+L4)*(h/6);
end
plot(t,x,'b',t,y,'r')
xlabel('time')
ylabel('x in blue and y in red')
figure
plot(x,y,'g')
title('2D figure(RK4)')
xlabel('X')
ylabel('Y')