%%% Note: This uses symbolic Matlab commands for the problems you are supposed to
%         do with pencil and paper. But if I ask you to write Matlab code in a problem
%         you should *NOT* use any symbolic commands.

function f460

disp('%%%%%%%%% Problem 1 %%%%%%%%%%')

t = sym('[-2; -1; 0; 1]')
y = sym('[ 2;  1; 2; 2]')
A = [abs(t), t.^2]
B = A'*A
b = A'*y
c = B\b

disp('%%%%%%%%% Problem 2 %%%%%%%%%%')

syms x1 x2 x3 real
x = [x1; x2; x3]
f = [4*x1 + x2^2 - 1; 4*x2 + x3^2 - 1; 4*x3 - x1^2 - 1]
Df = jacobian(f,x)

x0 = sym('[1;1;1]')
b = subs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)})
A = subs(Df,{x1,x2,x3},{x0(1),x0(2),x0(3)})
d = A\-b
xnew = x0 + d

disp('%%%%%%%%% Problem 3 %%%%%%%%%%')

%% see definition of function F at bottom

h = sym('1/2');
y0 = sym([1;2]); t0 = sym('1')

%% Euler
s = F(t0,y0);
y1 = y0 + h*s;
t1 = t0 + 1/2;
s = F(t1,y1);
y2 = y1 + h*s

%% RK2

s1 = F(t0,y0)
t1 = t0 + 1/2;
s2 = F(t1,y0+h*s1)
y1RK = y0 + h*(s1 + s2)/2

%%
format compact; format long g
for v=-2:.5:-1
    disp([v,solT(v)])
end

fzero(@solT,[-2,-1])

disp('%%%%%%%%% Problem 4 %%%%%%%%%%')

g = inline('1/x','x');
a = sym(1); b = sym(2);
h = (b-a)/2;
Q = h*(g(a)/2 + g(a+h) + g(b)/2)

%%% error formula: |Q-I| <= (1/12)*max|f''|*(b-a)^3/N^2
%  f = 1/x,   f' = -1/x^2,  f'' = 2/x^3,  max|f''(t)| = f''(1) = 1  for t in [1,2]
%  |Q-I| <= (1/12)*2*1^3/N^2 = 1/(6*N^2)
% we want 1/(6*N^2) <= 10^-15   <=>   N >= sqrt(10^15/6)
n = sqrt(10^15/6)
N = ceil(n)

%%%%%%%%%%%%% function F for problem 3 %%%%%%%%%
function yp = F(t,y)
yp = [y(2); t - y(1)*y(2)];

function z = solT(v)
[ts,ys] = ode45(@F,[0,2],[1;v]);
z = ys(end,1);


