runge-kutta ode45

I want to solve a two EDO system with the runge kutta method, I wrote this code but i don't get anything, someone know why this is wrong? if i use command ode45 is better

function [] = runge_kutta_sis(f,g,x0,y0,a,b,h)

t = a:h:b; n = length (t);

x = (x0); y = (y0);

for i=1:n-1

k1 = h*f(x(i),y(i),t(i));

l1 = h*g(x(i),y(i),t(i));

k2 = h*f(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);

l2 = h*g(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);

k3 = h*f(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);

l3 = h*g(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);

k4 = h*f(x(i)+k3,y(i)+l3,t(i)+h);

l4 = h*g(x(i)+k3,y(i)+l3,t(i)+h);

x (i+1) = x(i) + (1/6)*(k1+2*k2+2*k3+k4);

y (i+1) = y(i) + (1/6)*(l1+2*l2+2*l3+l4);

end

plot (t,x,t,y)

I write this in the command window but I can't see any figure

>> f = @(x,y,t) 10-0.1*x*y;

>> g = @(x,y,t) 20-0.2*x*y;

>> runge_kutta_sis(f,g,500,600,0,20,0.5)

Answer 1

[KALYAN ACHARJYA]

Minor modification: Function Code (Save in the same working directory, filename runge_kutta_sis.m)

function runge_kutta_sis(f,g,x0,y0,a,b,h)
t=a:h:b;
n=length(t);
x=zeros(1,length(n));
y=zeros(1,length(n));
x(1)=x0;
y(1)=y0;
for i=1:n-1
    k1 = h*f(x(i),y(i),t(i));
    l1 = h*g(x(i),y(i),t(i));
    k2 = h*f(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);
    l2 = h*g(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);
    k3 = h*f(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);
    l3 = h*g(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);
    k4 = h*f(x(i)+k3,y(i)+l3,t(i)+h);
    l4 = h*g(x(i)+k3,y(i)+l3,t(i)+h);
    x(i+1)=x(i)+(1/6)*(k1+2*k2+2*k3+k4);
    y(i+1)=y(i)+(1/6)*(l1+2*l2+2*l3+l4);
end
plot(t,x,t,y);
grid on;

Main Script:

f = @(x,y,t) 10-0.1*x*y;
g = @(x,y,t) 20-0.2*x*y;
runge_kutta_sis(f,g,500,600,0,20,0.5)

Output:

Hope it Helps!

[Alvaro Mª Zumalacarregui Delgado]

it helps thanks! if I want to do it with the command ode45 how can I write it?

[KALYAN ACHARJYA]

https://in.mathworks.com/help/matlab/ref/ode45.html

Answer 2

[Meysam Mahooti]

https://www.mathworks.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78?s_tid=srchtitle