function y=rk_hw5(tmin, tmax, nsteps, y0) % % User supplies initial time (tmin), final time (tmax), % number of time steps (nsteps), and vector of initial % values (y0). (y0 must be a COLUMN vector.) % This code uses the fourth-order Runge-Kutta method % to solve the initial value problem % % y'=f(t,y), y(tmin) = y0 % % The code returns an array y of size (d,nsteps+1), where d % is the number of unknowns (the length of y0), containing the % values of y at each timestep (including the initial time). % % Load initial values. % d = length(y0); y = zeros(d,nsteps+1); y(:,1) = y0; yk = y0; % % Loop over time steps. % h = (tmax-tmin)/nsteps; for k=1:nsteps, tk = tmin + k*h; F1 = f(tk, yk); F2 = f(tk+h/2, yk+(h/2)*F1); F3 = f(tk+h/2, yk+(h/2)*F2); F4 = f(tk+h, yk+h*F3); yk = yk + (h/6)*(F1 + 2*F2 + 2*F3 + F4); y(:,k+1) = yk; end; % % Function defining y'. % function fty = f(t,y) fty(1,1) = y(3); fty(2,1) = y(4); fty(3,1) = -y(1)/(y(1)^2 + y(2)^2)^(3/2); fty(4,1) = -y(2)/(y(1)^2 + y(2)^2)^(3/2); return;