% Suppose for the same problem from lecture 16, of y'' = -y, 
% with boundary conditions y(0) = -1 and y(1) = 2

% In the general ODE from lecture 17 that would imply:
% p(t) = 0, q(t) = -1 and r(t) = 0.

% Resulting in 
% A(t) = 1
% B(t) = -[2 + q(t)dt^2]
% C(t) = 1
% rhs = 0;

% Suppose we want this for N times tn, then we
% construct a matrix of size N-2:

N = 10;        % number of total points
n = N-2;       % number of interior points
dt = 1/(n+1);  % time step size

% Construct the matrix A using the function diag, or spdiags (for sparse
% computations
diagonal_vec = -(2+(dt^2)*(-1))*ones(n,1);
diagonal_plus_one_vec = ones(n-1,1);
diagonal_minus_one_vec = ones(n-1,1);
  
A = diag(diagonal_vec,0) + diag(diagonal_plus_one_vec,1) + diag(diagonal_minus_one_vec,-1);
b = zeros(n,1);
b(1) = 0 - 1*(-1);
b(end) = 0 - 1*2;

% Solve the system of equations
x = A\b;

% Plot solution including initial point and end points: (0,-1) and (1,2)
t = (0:dt:1)';
sol = zeros(N,1);
sol(1) = -1;
sol(end) = 2;
sol(2:N-1) = x(1:n);
plot(t,sol,'*-'); hold on;


% Plot exact solution (from analytic solution of ODE)
C1 = (2*exp(i) + 1)/(exp(2*i)-1);
C2 = -1 - C1;
plot(t,real(C1*exp(i*t) + C2*exp(-i*t)),'k^-');