% 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
                                % update x directly.
end
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. I inline the modified version and remove any comments to keep
% mine in focus.

% f = inline('y-t^2+1','t','y'); % We have to replace this one with a function
                                 % describing our system of equations.
t0 = 0;
tend = 10;                       % Change integration interval.
y0 = [1; 2];                     % New initial value - this time a vector.

Nstep = 100;                     % New number of steps.
h = (tend-t0)/Nstep;

t = t0;
y = y0;

yn = y0;
tn = t0;
for n = 1:Nstep
    yn = yn+h*ourfunction(tn,yn); % Call our new function.
    tn = tn+h;
    y = [y,yn];
    t = [t,tn];
end

plot(t,y,'-o')

% 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);

