% Task 1

% First we input the matrix and the vector.
A = [2, 1, 1; 2, 3, 2; 1, 1, 2]
b = [7; 7; 6]

% To solve the system in MATLAB we simply use the '\' operator.
x = A \ b;

% (Visual) check that the answer is correct, we simply print Ax.
disp('Should be equal to b:')
A * x

% We obtain the eigenvectors and eigenvalues using the command 'eig'.

[V, D] = eig(A);

disp('Eigenvectors of A:')
V
disp('Eigenvalues of A:')
diag(D)


% Task 2

% We inline this here for clarity of presentation. I really recommend that you
% use separate files per function/script.

disp('Approximation of sqrt(3) using 5 iterations of Newton''s method:')

x = 2; 							    % Initial guess.
for i = 1 : 5						% We do 5 iterations.
    x = x - (x^2 - 3) / (2 * x);	% Since we do not need the previous values, we
end                                 % update x directly.

x

% We then want to add a stopping criterion. We want to stop when the difference
% between two iterates is less than 10^(-6), i.e. when we have at least 6 correct
% digits after the comma in the answer. The main thing we have to do is exchange
% the 'for' with a 'while'.

disp('Approximation of sqrt(3) using an error tolerance of 10^(-6):')

x   = 2;			    			% Initial guess.
err = 1;			    			% Initial error > tolerance.
j	= 0;		    				% Number of iterations.
while err > 10^(-6) && j < 100 	    % We want to iterate until the error is smaller
                                    % than 10^(-6). In addition we put a maximum of
                                    % 100 iterations to avoid the program from 
                                    % freezing up if something nasty happens.

    y   = x - (x^2 - 3) / (2 * x);	% Since we need the previous value of x, we use
                                    % a temporary variable to store the new iterate.

    err = abs(y - x);    			% Calculate the iterate change.
    x   = y;		       			% Throw away the old value.
    j	= j + 1;	    			% Increase the number of iterations by 1.
end
x
disp('Number of iterations:')
j


% Task 3

% You were given a file 'euler.m' on the home page. I inline it here for clarity
% of presentation.

% function [tn, yn] = euler(f, t0, tend, y0, h)
% %
% % Solve the ODE y' = f(t, y) with Eulers method. 
% % Input: f       : Function handle to the rhs. of the ODE. 
% %        t0, tend: The integration interval. 
% %        y0      : Initial value.   
% %        h       : Stepsize.
% % Output: tn     :
% %         yn     : The solution at the points tn.
% 
% t = t0;    % Current solution
% y = y0;
% 
% tn = t0;   % Store the solution
% yn = y0;
% 
% while t < tend
%     if t + h > tend;
%         h = tend - t;
%     end
%     
%     % The Euler step
%     y = y + h * feval(f, t, y);
%     t = t + h;
%     
%     % Store the solution
%     tn = [tn, t];
%     yn = [yn, y];
% end

t0 = 0;                                             % Start time
tend = 10;                                          % End time
y0 = [1; 2];                                        % Initial value
h = 0.1;                                            % Step size

[tn, yn] = euler(@ourfunction, t0, tend, y0, h);    % Solve using euler.m

plot(tn, yn, '-o')                                  % Plot the result

% ourfunction has to be saved in a separate file, ourfunction.m, inlined here
% for clarity.

% function yp = ourfunction(t, y)     % A function named ourfunction taking
%                                     % two parameters (t, y) and returning one,
%                                     % yp.
%     yp    = zeros(2, 1);             % To assure that we get a column vector.
%     yp(1) = y(1) * (1 - y(2));
%     yp(2) = y(2) * (y(1) - 2);
% end

