


a=0;b=pi;c=1;
sol=@(x) sin(x);
f=@(x) (1+c)*sin(x);
N=10;
M=5;
err=zeros(M,2);
%equidistant mesh
for k=1:M,
  x=linspace(a,b,N+1)';
  [A,B,h]=FEM1Dmats(x);
  %Since we have dirichlet bnd cond we cat extract submatrices
  iN=2:N;
  iA=A(iN,iN);iB=B(iN,iN);
  ix=x(iN);
  %FDM and FEM 1D RHS (same)
  F1=F2=f(ix);
  u2=(iA+c*iB)\(h*F2); %solves (iA+c*iB)*u2=F2
  u1=(iA+c*(h)*eye(size(iB)))\(h*F1);
  plot(ix,u1,ix,u2,x,sol(x));pause(1);
  er1=u1-sol(ix);
  er2=u2-sol(ix);
  err(k,:)=[norm(er1,"inf"),norm(er2,"inf")];
  N*=2;
endfor
err(1:M-1,:)./err(2:end,:)


a=0;b=pi;c=1;
sol=@(x) sin(x);
f=@(x) (1+c)*sin(x);
N=10;
M=6;
err=zeros(M,2);
%non-equidistant  mesh
for k=1:M,
  x1=linspace(a,b,N+1)';
%  x2=pi*sin(linspace(0,1,N+1))/sin(1);
  x2=x1;
  x2(3)=0.5*(x2(2)+x2(3));
  [A1,B1,h1]=FEM1Dmats(x1);
  [A2,B2,h2,hh2]=FEM1Dmats(x2);
  %Since we have dirichlet bnd cond we cat extract submatrices
  iN=2:N;
  iA1=A1(iN,iN);
  iA2=A2(iN,iN);iB2=B2(iN,iN);
  ix1=x1(iN);
  ix2=x2(iN);
  %FDM and FEM 1D RHS (same)
  F1=h1*f(ix1);
  F2=f(ix2).*(hh2(1:N-1)+hh2(2:N))*0.5;
  u2=(iA2+c*iB2)\(F2); %solves (iA+c*iB)*u2=F2
  u1=(iA1+c*(h1)*eye(size(iA1)))\(F1);
  %plot(ix1,u1,ix2,u2,x,sol(x));pause(1);
  er1=u1-sol(ix1);
  er2=u2-sol(ix2);
  err(k,:)=[norm(er1,"inf"),norm(er2,"inf")];
  N*=2;
endfor
err
err(1:M-1,:)./err(2:end,:)



a=0;b=pi;c=1;
sol=@(x) sin(x);
f=@(x) (1+c)*sin(x);
NN=5;
M=5;
err=zeros(M,4);
%non-equidistant  mesh on pi*[0.5,1] mesh dobled
for k=1:M,
  x=linspace(a,0.5*(a+b),NN+1);
  x2=linspace(0.5*(a+b),b,2*NN+1);
  x2=[x,x2(2:end)];x2=x2';
  x1=linspace(a,b,3*NN+1)';
  [A1,B1,h1,hh1]=FEM1Dmats(x1);
  [A2,B2,h2,hh2]=FEM1Dmats(x2);
  N=3*NN;
  %Since we have dirichlet bnd cond we cat extract submatrices
  iN=2:N;
  iA1=A1(iN,iN);
  iA2=A2(iN,iN);iB2=B2(iN,iN);
  ix1=x1(iN);
  ix2=x2(iN);
  %FDM and FEM 1D RHS (same)
  F1=h1*f(ix1);
  F2=f(ix2).*(hh2(1:N-1)+hh2(2:N))*0.5;
  u2=(iA2+c*iB2)\(F2); %solves (iA+c*iB)*u2=F2
  u1=(iA1+c*(h1)*eye(size(iA1)))\(F1);
  %plot(ix1,u1,ix2,u2,x,sol(x));pause(1);
  er1=u1-sol(ix1);
  er2=u2-sol(ix2);
  err(k,:)=[3*N,h1,norm(er1,"inf"),norm(er2,"inf")];
  NN*=2;
endfor
err
err(1:M-1,:)./err(2:end,:)
