% Script to solve HW3 problems 1b-c

% Backward-Euler and BDF2 on 
%  dA/dt = - lambda*A + 1 - cos(t)
%  dB/dt =   lambda*A - B
%  A(0) = B(0) = 0. Want B(T) at T = pi.
%  If we write phi = [A; B], this pair of ODEs has the form
%  dphi/dt = M*phi + s
%   where M = [-lambda 0; lambda -1] and s = [1 - cos(t); 0]
  lambda = 1000;
  T = pi;
  M = [-lambda  0; lambda -1];
  
  pmin = 2;
  pmax = 7;
  np = pmax - pmin + 1;

  for ip = 1:np
    N = 2^(pmin+ip-1);
    dt = T/N;

    % Backward Euler
    %  (1 - M*dt)phi(n+1) = phi(n) + dt*s(n+1)
    
    phi = [0; 0];
    for it = 1:N
      tnp1 = it*dt;
      s = [1 - cos(tnp1); 0];
      phi = (eye(2) - M*dt)\(phi+dt*s);
    end
    B1(ip) = phi(2); % B at T=pi

    % BDF2

    % Starting Backward-Euler step
    lamdt = lambda*dt;
    phim1 = [0; 0];
    tnp1 = dt;
    s = [1 - cos(tnp1); 0];
    phi = (eye(2) - M*dt)\(phim1+dt*s);
    
    % BDF2 steps:
    % (3 - 2*M*dt) phi(n+1) = 4*phi(n) - phi(n-1) + 2*dt*S(n+1)

    for it = 2:N
      tnp1 = it*dt;
      s = [1 - cos(tnp1); 0];
      phip1 = (3*eye(2) - 2*M*dt)\(4*phi - phim1 + 2*dt*s);
      phim1 = phi;
      phi = phip1;
    end
    B2(ip) = phi(2); % B at T=pi
  
  end
  incr1 = abs(B1(2:np)-B1(1:np-1));
  incr2 = abs(B2(2:np)-B2(1:np-1));

  dxrange = T./2.^(pmin+1:pmax);
  loglog(dxrange, incr1, 'bx',dxrange, incr2, 'r+',...
       dxrange,incr1(np-1)*dxrange/dxrange(np-1), 'b-',...
       dxrange,incr2(np-1)*(dxrange/dxrange(np-1)).^2, 'r-')
  xlabel('\Delta t')
  ylabel('\epsilon')
  legend('B-E','BDF2','\Delta t^1 fit line','\Delta t^2 fit line',4)