function RB = Solve_RK4(RB)

K = length(RB);

dt = get(RB{1},'dt');

state = cell(K,1);
k1 = cell(K,1);
k2 = cell(K,1);
k3 = cell(K,1);
k4 = cell(K,1);
R = cell(K,1);
for k = 1:K
  state{k} = zeros(12,1);
  k1{k} = zeros(12,1);
  k2{k} = zeros(12,1);
  k3{k} = zeros(12,1);
  k4{k} = zeros(12,1);
  R{k} = zeros(12,1);
end

N = get(RB{1},'N');
for n = 1:N-1
  disp(['Timestep = ' num2str(n)]); 

  for k = 1:K
   state{k}(1:6) = get_pos(RB{k},n);
   state{k}(7:12) = get_vel(RB{k},n);
   RB{k} = set(RB{k},'sv',state{k});
  end
  for k = 1:K
   [rhs F] = RHS(RB{k},RB,n);
   k1{k} = dt*rhs;
   for i = 1:RB{k}.n_Forces
     RB{k} = set_force(RB{k},F{i},i,n);  %Store force results at this timestep.
   end
  end
  
  for k = 1:K
   RB{k} = set(RB{k},'sv',state{k}+k1{k}/2);
  end
  for k = 1:K
   [rhs F] = RHS(RB{k},RB,n+.5);    
   k2{k} = dt*rhs;
   
  end

  for k = 1:K
   RB{k} = set(RB{k},'sv',state{k}+k2{k}/2);
  end
  for k = 1:K
   [rhs F] = RHS(RB{k},RB,n+.5);    
   k3{k} = dt*rhs;
  end

  for k = 1:K
   RB{k} = set(RB{k},'sv',state{k}+k3{k});
  end
  for k = 1:K
   [rhs F] = RHS(RB{k},RB,n+1);       
   k4{k} = dt*rhs;
  end

  for k = 1:K
   R{k} = k1{k}/6+k2{k}/3+k3{k}/3+k4{k}/6;
   state{k} = state{k} + R{k};
  end

  for k = 1:K
   RB{k} = set_pos(RB{k},state{k}(1:6),n+1);
   RB{k} = set_vel(RB{k},state{k}(7:12),n+1);
   RB{k} = set_accel(RB{k},R{k}(7:12)/dt,n+1);
  end

  
end

%Compute and store force for last timestep.
for k = 1:K
   state{k}(1:6) = get_pos(RB{k},N);
   state{k}(7:12) = get_vel(RB{k},N);
   RB{k} = set(RB{k},'sv',state{k});
end
for k = 1:K
 [rhs F] = RHS(RB{k},RB,N);
 for i = 1:RB{k}.n_Forces
   RB{k} = set_force(RB{k},F{i},i,N);  %Store force results at this timestep.
 end
end
