Blatt6_erste

function r_mat = test2(y0,t_vec)
  h=t_vec(2)-t_vec(1);
  r_mat= zeros(length(y0),length(t_vec));
  r_mat(:,1)=y0;
  for i=2:length(t_vec)
   k2=r_mat(:,i-1)+h*2/3*f(t_vec(i-1),r_mat(:,i-1));
   k3=r_mat(:,i-1)+h*2/3*f(t_vec(i-1)+2/3*h,k2);
   r_mat(:,i)=r_mat(:,i-1) +0.25*h*f(t_vec(i-1),r_mat(:,i-1))+0.375*h*f(t_vec(i-1)+2/3*h,k2)+0.375*h*f(t_vec(i-1)+2/3*h,k3);
  endfor


endfunction



function z = f(t,y)
  m1=1;
  m2=10;
  m3=10;

  a1=m2/(norm(y(4:6)-y(1:3))^3);
  b1=m3/(norm(y(7:9)-y(1:3))^3);

  a2=m2/(norm(y(1:3)-y(4:6))^3);
  b2=m3/(norm(y(7:9)-y(4:6))^3);

  a3=m2/(norm(y(1:3)-y(7:9))^3);
  b3=m3/(norm(y(4:6)-y(7:9))^3);

  A=zeros(length(y),length(y));
  A(1:9,10:18)=eye(9);
  A(10:12,1:3)=-(a1+b1)*eye(3);
  A(10:12,4:6)=(a1)*eye(3);
  A(10:12,7:9)=(b1)*eye(3);

  A(13:15,4:6)=-(a2+b2)*eye(3);
  A(13:15,1:3)=(a2)*eye(3);
  A(13:15,7:9)=(b2)*eye(3);

  A(16:18,7:9)=-(a3+b3)*eye(3);
  A(16:18,1:3)=(a3)*eye(3);
  A(16:18,4:6)=(b3)*eye(3);
  A;

  z= A*y;

endfunction

%{
Startvektor= [1 0 0 0 1 0 -1 0 0 0 -1 0 0 0 0 0 0 1].';
test2(Startvektor,[0:0.01:1])
plot3(ans(1,:),ans(2,:),ans(3,:),'color','r',ans(4,:),ans(5,:),ans(6,:),'color','g',ans(7,:),ans(8,:),ans(9,:),'color','b')
%}