function err=refined_err(x,y)
% err=refined_err(x,y)
% compute error using a refined mesh
% x is the nodal mesh, including endpoints
% Y is the calculated solution at the nodes
% use phi.m to interpolate values of solution at intermediate nodes
% Objective is to show that convergence is really O(h^3) in Exercise 7
% intended to be run with supporting files from Lab 4

% mms 3/8/17

% N is number of nodes, including endpoints
N=numel(x);
if mod(N,2)==0
  error('N MUST be odd for quadratic elements!')
end
h=x(2)-x(1);

% compute nodes for refined mesh
% with midpoints, total of N+N-1 nodes
x=reshape(x,N,1); % make x a column vector
xmid=.5*(x(2:N)+x(1:N-1));
xref=zeros(2*N-1,1);
xref(1:2:2*N-1)=x;
xref(2:2:2*N-2)=xmid;

% not assuming phi works with vector inputs
ymid(1)=y(1)*phi(0,h,xmid(1))+y(2)*phi(1,h,xmid(1))+y(3)*phi(2,h,xmid(1));
ymid(2)=y(1)*phi(0,h,xmid(2))+y(2)*phi(1,h,xmid(2))+y(3)*phi(2,h,xmid(2)) ...
       +y(4)*phi(3,h,xmid(2));
for kphi=2:N-3  % k indexes coarse mesh nodes and phis
  ymid(kphi)=y(kphi-1)*phi(kphi-2,h,xmid(kphi)) ...
      +y(kphi)*phi(kphi-1,h,xmid(kphi)) ...
      +y(kphi+1)*phi(kphi,h,xmid(kphi)) ...
      +y(kphi+2)*phi(kphi+1,h,xmid(kphi)) ...
      +y(kphi+3)*phi(kphi+2,h,xmid(kphi));
end
kphi=N-2;
  ymid(kphi)=y(kphi-1)*phi(kphi-2,h,xmid(kphi)) ...
      +y(kphi)*phi(kphi-1,h,xmid(kphi)) ...
      +y(kphi+1)*phi(kphi,h,xmid(kphi)) ...
      +y(kphi+2)*phi(kphi+1,h,xmid(kphi));
kphi=N-1;
  ymid(kphi)= ...
      +y(kphi-1)*phi(kphi-2,h,xmid(kphi)) ...
      +y(kphi)*phi(kphi-1,h,xmid(kphi)) ...
      +y(kphi+1)*phi(kphi,h,xmid(kphi));

yref=zeros(2*N-1,1);
yref(1:2:2*N-1)=y;
yref(2:2:2*N-2)=ymid;

err=max(abs(yref-exact7(xref)));

%{
For sizes [7,15,31,61,121],
err =    1.2321e-04   1.5377e-05   1.9161e-06   2.6291e-07   3.4465e-08
rats =    8.0125    8.0251    7.2881    7.6283
%}
