%Subject: Problem Set #16
%% Problem 1
clear all  clc
A = [1, 3, -5, 9, -1; 2, -2, -3, -7, 2; 1, -5, 2, -16, 3]; B = [5; 5; 0];
A_a = [A B];
sol = rref(A_a)
sol_symb = sym(sol)
disp('There are infinite many roots')
%% Problem 3
clc, format longG
x_int = [0,1];
y0 = [1;6];
sol = ode45(@syst_prob, x_int, y0);
x = 0:0.01:1;
y = deval(sol, x);
plot(x,y(1,:), 'r', 'LineWidth', 3)
legend('y(x)')
xlabel('x'), ylabel('y(x)')
figure
plot(x, y(1,:), 'r', x,y(2,:), 'g', 'LineWidth', 3)
legend('First order derivative','Second order derivative')
xlabel('x'), ylabel('yprime')

function dydx = syst_prob(x,y)
    dydx = zeros(2,1);
    dydx(1) = y(2);
    dydx(2) = 4*x*y(2)-6*y(1)/1+x^2;
end

%% Problem 4
clc, clear all,close all
format compact, format longG
izt=@(a,x)(exp(a*(x.^1/2)))/4 + sin(x);
for a=1:10
 fun=@(x)izt(a,x);
 int_value(a)=integral(fun,0,1) - 2;
end
a_val=1:10;
plot(a_val,int_value,'r','LineWidth',3)
grid
xlabel('a'),ylabel('f(a)')
a_p=[5.8 5.9]; iter=6;
for i=1:2
 fun=@(x)izt(a_p(i),x); f(i)=integral(fun,0,1) - 2;
end
for k=1:iter
 k1=k+1; i=k+2;
 a_p(i)=a_p(k1)-f(k1)*(a_p(k1)-a_p(k))/(f(k1)-f(k));
 fun=@(x)izt(a_p(i),x); f(i)=integral(fun,0,1)-2;
 M_pr(k,1)=a_p(i);M_pr(k,2)=f(i);
end
disp('approximation to the root:'), disp(a_p(1:2))
disp('               a                    f(x)'), disp(M_pr)
disp('The value of the root is: 5.9107')
%% Problem 6
clear all, clc, format compact
A = [15, 10, 9; 10, 12, 10; 9, 10, 12]; B = [5; 21; 2];
xapp = zeros(3,1);
iter = 0;
itermax = 60;
k = 1;
while iter < itermax
    k = k+1; iter = iter+1;
    for i = 1:3
        rez = 0;
        for j = 1:3
            if j~=i
                rez=rez+xapp(j, k-1)*A(i,j);
            end
        end
        xapp(i,k) = (B(i,1)-rez)/A(i,i);
    end
end
x_approx = xapp(:, k)
x_sol = linsolve(A,B)

%% Problem 5
clear all, clc
A = [7, 7, 9; 7, 16, -7; 9, -7, 107]; B = [19;37;4];
ni = 1;
for i = 1:3
    if det(A(1:i, 1:i))>0
    else ni = 2; break
    end
end
if ni == 2
    disp('the coefficient matrix is not positive definite'), return
end
check=isequal(A,A');
if check == 0
    disp('the coefficient matrix is not symmetric'), return
end
disp('A  is symmetric and positive definite ')

k_iter = 0; epsi = 10^(-3); itermax = 300; x = [0 2 0];
r = A*x-B; norm_r = norm(r);
while norm_r > epsi & k_iter < itermax
    k_iter=k_iter+1;
    tau=((A*r)'*r)/norm(A*r)^2;
    x=x-(tau*r')';r=A*x-B; norm_r=norm(r);
end
k_iter,tau,x,norm_r
x_sol=linsolve(A,B)
disp('Answer:'),

%% Problem 2
clear all, close all, clc
graph1 = ezplot('x1.^3+x2.^3-5',[-1 2 0 3]) %assume that x=x1 & y=x2
hold on
graph2 = ezplot('sin(x2+0.8)-x1-0.7',[-1 2 0 3])
set(graph1,'Color','r','LineWidth',3)
set(graph2,'Color','b','LineWidth',3)
hold off
grid

syms x1 x2
xapp=[-0.8 1.70]; xapp_pr=xapp; xpr=[x1 x2];
fun=[x1.^3+x2.^3-5,sin(x2+0.8)-x1-0.7]; % f =(f1 f2)
fun_pr=[diff(fun(1),xpr(1)) diff(fun(1),xpr(2))
 diff(fun(2),xpr(1)) diff(fun(2),xpr(2))]% partial derivatives

epsi=10^(-5); k=1; % number of iterations
norm_sol=1;
while norm_sol > epsi
 for i=1:2
 B(i,1)=-double(subs(fun(i),xpr,xapp));
 for j=1:2
 A(i,j)=double(subs(fun_pr(i,j),xpr,xapp));
 end
 end
 xdelta=A\B; xapp=xapp+xdelta';
 c=double(subs(fun,xpr,xapp)); norm_sol=norm(c);
 M_pr(k,1)=xapp(1); M_pr(k,2)=xapp(2); M_pr(k,3)=norm_sol; k=k+1;
end

disp('approximation to the root:'), disp(xapp_pr(1:2))
disp('      x1      x2       error norm')
disp(M_pr)