%% Numerik D-MAVT, RW (Nipp), SS05
%
% This is the M-file solution to the twelveth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
%
% Then to run the file type "s12" in the command window

function s11
% Aufgabe 1
disp(['Aufgabe 1']);
A=[-100 1;0 -1/10]; 
y0=[1;999/10];
l=-1/10;
I=eye(size(A));
% a), c)
h=0.1;
yE = y0; yT = y0;
eE = 0.; eT = 0.;	% error at t=0
ME = I+h*A; MT = inv(I-h/2*A)*(I+h/2*A);
for t=h:h:1.8
  yE = ME*yE;
  yT = MT*yT;
  y  = y0*exp(l*t);
  eE = [eE norm(yE-y)];
  eT = [eT norm(yT-y)];
end
figure;
subplot(2,1,1);
plot(eE); title('Euler explizit mit h=0.1');
ylabel('error in 2 norm'); xlabel('n');
subplot(2,1,2);
plot(eT); 
title('Trapez mit h=0.1');
ylabel('error in 2 norm'); xlabel('n');

% c)  
hh=[0.015 0.02 0.025];
e = [];
for h=hh
  ME = I+h*A;
  yE = y0;
  e = [e 0];
  for n=1:80
    yE = ME*yE; 
    y  = y0*exp(l*n*h);
    e = [e norm(yE-y)];
  end
end
le=length(e);
figure;
plot(e(1:le/3),'--'); hold on; 
plot(e(le/3+1:2*le/3),'-.'); hold on;
ylabel('error in 2 norm'); xlabel('n');
title('Euler explizit');
legend('h<h_{krit}','h=h_{krit}',2);


%%%%%% AUFGABE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

disp(['Aufgabe 3 working...this may take some time:-)'])

a=7; % parameter for ellipse 
e=0.9;
y1=a*(e+1);
y2=0;
y3=0;
y4=0.2;
y0=[y1;y2;y3;y4];tf=10; 	% initial condition, final time
y0_halfstep=[y1;y2;y3;y4];


%h=0.03;
h=0.05;


h_halfstep=h/2;
t=0;
  
  y = y0; orbit=y(1:2);
  y_halfstep = y0_halfstep; orbit_halfstep=y_halfstep(1:2);

  TOL=0.00001;

  for i=1:15000

        y_try = y + h*feval('kepler',h,y);
        
        y_halfstep_try = y +h/2*feval('kepler',h/2,y);
        y_halfstep_try = y_halfstep_try +h/2*feval('kepler',h/2,y_halfstep_try);
        
        l=abs(y_try(3)-y_halfstep_try(3));
        m=abs(y_try(4)-y_halfstep_try(4));
        n=abs(y_try(1)-y_halfstep_try(1));
        o=abs(y_try(2)-y_halfstep_try(2));
        
     if (l+m+n+o)>=TOL
          h =h*(TOL/(l+m+n+o))*0.8;
          y_try = y + h*feval('kepler',h,y);
          y_halfstep_try = y +h/2*feval('kepler',h/2,y);
          y_halfstep_try = y_halfstep_try +h/2*feval('kepler',h/2,y_halfstep_try);
          l=abs(y_try(3)-y_halfstep_try(3)); 
          m=abs(y_try(4)-y_halfstep_try(4));
          n=abs(y_try(1)-y_halfstep_try(1));
          o=abs(y_try(2)-y_halfstep_try(2));
     end
     step(i)=h; % store step h for plot
      t=t+h;
      time(i)=t; % store time for plot
      h=h*(TOL/(l+m+n+o))*0.8;
      y = y_try;
      orbit = [orbit y(1:2)];
    %disp(time(i));
    %disp(i);

  end
  figure;
  plot(time,step);
  xlabel('t');ylabel('h(t)');
  title(['h(t) fuer Euler explizit mit Schrittweitensteuerung']);
  disp(t);
  figure;  
  plot(orbit(1,:),orbit(2,:),'.:');
  xlabel('y_1(t)');ylabel('y_2(t)');
  title(['Aufgabe 3: Kepelr motion with Euler expicit']);
  
%--------------------------------------------------------------------------
function f=kepler(t,y)

d=(y(1)*y(1)+y(2)*y(2))^(1.5);
f=[y(3);y(4);-y(1)/d;-y(2)/d];
%------------------------------------------------------------------------







