function [t,y] = RungeKutta4(fun, tspan, y0, At)
global MBody
n = length(MBody.Ind);
global joy
t0 = tspan(1);
tfin = tspan(2);
t(1) = t0;
y(:,1) = y0;
i = 2;
tt = t0;
yy = y0;
while tt<=tfin,
k1 = feval(fun,tt,yy);
k2 = feval(fun,tt+At,yy+(At)*k1);
yy = yy + (At/2)*(k1+k2);
tt = tt+At;
t(i) = tt;
y(:,i) = yy;
q = y(1:n);
i = i+1;
if button(joy, 1) == 1
break
end
end
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
