clear
clf

error=zeros(1,8);
f=@(x)(-6*pi*cos(3*pi*x) + 9*pi*pi*x.*sin(3*pi*x));

% loop over number of nodes
for k=3:8
  % nodes and mesh width  
  N = 2^k;
  h = 1/N;
  
  % prepare sparse matrix
  % first row
  A=(1/h^2)*spdiags([-ones(N-1,1) 2*ones(N-1,1) -ones(N-1,1)],-1:1,N-1,N-1);
  
  x=h:h:1-h;

  % right hand side
  b=f(x);
    
  % compute solution
  % direct solver
    u = A\b';
  % Jacobi
    res = 1;
    ite = 0;
    u = zeros(N-1,1);
    J=diag(A);
    while (res > 1e-10)
        res_vec = b' - A*u;
        u = u + res_vec./J;
        res = norm(res_vec);
        ite = ite + 1;
    end
   
    
  % plot solution
  plot(x,u,'-bd');
  
  % compute error
  u_ana = transpose(x.* sin(3*pi*x));
  error(k) = norm(u_ana-u)/sqrt(N-1);
  s = sprintf('nodes %d error %f order %f iterations %d',N,error(k),log(error(k-1)/error(k))/log(2),ite);
  disp(s)
end