% Kod för tentan 16/1 2015

%%
% Del 2, uppgift 1 b)

N = 100;
h = 1/(N+1);

k = @(x) 1+x;
kp = @(x) 1;
f = @(x) x;

A = zeros(N,N);
F = zeros(N,1);

for n=2:N-1
  x = n*h;

  A(n,n-1) = h*kp(x)/2 - k(x);
  A(n,n) = 2*k(x);
  A(n,n+1) = -h*kp(x)/2 - k(x);

  F(n) = h^2 * f(x);
end

A(1,1) = 2*k(h);
A(1,2) = -h*kp(h)/2 - k(h);

A(N,N-1) = h*kp(N*h)/2 - k(N*h);
A(N,N) = 2*k(N*h);

u=A\F;
figure(1)
plot(u)

%%
% Del 2, uppgift 2 a)

N = 1000000;
x = rand(N,1);
I = (1/N)*sum(exp(-x.^2))

%%
% Del 2, uppgift 2 b)

N = 1000000;
x = rand(N,3);
tmp = exp(x.^2);
tmp = 1./sum(tmp,2);
I = (1/N)*sum(tmp)

%%
% Del 2, uppgift 3 b)

b = 10;
N = 10000;
dt = b/N;

g = 9.8;
L = 1;

T = 0:dt:b;
X = zeros(2, N+1);
X(:,1) = [0;1];

for n=1:N
  X(:,n+1) = X(:,n) + dt*[X(2,n); -g/L*sin(X(1,n))];
end
figure(2)
plot(T,X(1,:))

%%
% Del 2, uppgift 3 c)

b = 10;
N = 10000;
dt = b/N;

g = 9.8;
L = 1;

T = 0:dt:b;
X = zeros(2, N+1);
X(:,1) = [0;1];

for n=1:N
  X(:,n+1) = X(:,n) + dt*[X(2,n); -g/L*sin(X(1,n))];
end
plot(T,X(1,:))

% Tillägg:
for n=1:N
  if X(1,n) < 0 && X(1,n+1) > 0
    Tperiod = T(n)
    break
  end
end
