Animalitos

%Assign values
n = 2;
xi = 0;
xf = 2;
yi(1) = 4;
yi(2) = 6;
dx = 0.5;
xout = 0.5;

m = 1;
x = xi;
%set initial conditions
for i =1 : n
y(i) = yi(i);
end
%save output values
xp(m) = x;
for i = 1 : n
yp(m, i) = y(i);
end
while (1) %Print loop
xend = x+ xout;
if (xend > xf)
xend = xf; %Trim step if increment exceeds
end
h = dx;
m = m + 1;
%save output valuesxp(m) = x;
for i = 1 : n
yp(m, i) = y(i);
end
if (x >= xf), break, end
end

%Sub Integrator(x, y, h, n, xend)

while (1)%Calculation loop
if (xend - x < h)
h = xend - x; %Trim step if increment exceeds
end
%Call RK4Sys(x, y, h, n, ynew)
forj = 1 : n
y(j) = ynew(j);
end
if (x >= xend) , break, end
end
%End Sub

%Sub RK4Sys(x, y, h, n, ynew)

%Implement RK4 methodfor systems of ODEs
%Call Derivs(x, y, k1)
for j = 1 : n
ym(j) = y(j) + k1(j) * h / 2;
end
%Call Derivs(x + h / 2, ym, k2)
for j = 1: n
ym(j) = y(j) + k2(j) * h / 2;
end
%Call Derivs(x + h / 2, ym, k3)
for j = 1 : n
ye(j) = y(j) + k3(j) * h;
end
%Call Derivs(x +h, ye, k4)
for j = 1 : n
slope(j) = (k1(j) + 2 * (k2(j) + k3(j)) + k4(j)) / 6;
end
for j = 1 : n
ynew(j) = y(j) + slope(j) * h;
endx = x + h;
%End Sub

%Sub Derivs(x, y, dydx)
%Define ODE
dydx(1) = -0.5 * y(1);
dydx(2) = 4 - 0.3 * y(2) - 0.1 * y(1);
%End Sub