%% Sauer 6.3.1

h = 1/4;
t_end = 1;
n_steps = t_end/h;
% a

fA = @(y) [y(1)+y(2);-y(1)+y(2)];
fB = @(y) [-y(1)-y(2);y(1)-y(2)];
fC = @(y) [-y(2);y(1)];
fD = @(y) [y(1)+3*y(2);2*y(1)+2*y(2)];

yA = zeros(2,n_steps+1);
yB = zeros(2,n_steps+1);
yC = zeros(2,n_steps+1);
yD = zeros(2,n_steps+1);

yA(:,1) = [1;0];
yB(:,1) = [1;0];
yC(:,1) = [1;0];
yD(:,1) = [5;0];

yA_ana = [exp(1)*cos(1);-exp(1)*sin(1)];
yB_ana = [exp(-1)*cos(1);exp(-1)*sin(1)];
yC_ana = [cos(1);sin(1)];
yD_ana = [3*exp(-1)+2*exp(4);-2*exp(-1)+2*exp(4*1)];

for i = 1:n_steps
   yA(:,i+1) = yA(:,i) + h * fA(yA(:,i));
   yB(:,i+1) = yB(:,i) + h * fB(yB(:,i));
   yC(:,i+1) = yC(:,i) + h * fC(yC(:,i));
   yD(:,i+1) = yD(:,i) + h * fD(yD(:,i));
end

errorA = abs(yA(:,end)-yA_ana)
errorB = abs(yB(:,end)-yB_ana)
errorC = abs(yC(:,end)-yC_ana)
errorD = abs(yD(:,end)-yD_ana)
;