
function [u,x,B,F]=lap1mix17(f,a,b,al,be,N,mix=1,c=1)
%FDM solver for -u''+cu=f in (a,b)
%with u(a)=al u(b)=be (Dirichlet or mixed bnd cond)
%INPYT
%f- function handle t ooctave's function computing f(x)
%   (it should work for vectors i.e. f(z)=(f(z(k)) if z a vector
%a,b - ends of interval
%al,be - left/right bnd value
%N - mesh points no minus one
%mix - 0 Dirichlet bnd,  -1 at left bnd point i.e .a 
%      Neumann bnd; 1 at right bnd point N. bnd 
%c - constant coefficient (default =0)
%OUTPUT 
%u - FDM soltion on equidistant mesh
%x- mesh
%EXAMPLE
%f=sol=@sin;dsol=@cos;c=0;a=0;b=1;al=sol(a);be=sol(b);mix=0;
%N=20;[u,x]=lap1mix17(f,a,b,al,be,N,mix=0,c);er=u-sol(x);norm(er,'inf')
%N=20;[u,x]=lap1mix17(f,a,b,dsol(a),be,N,mix=-1,c);er=u-sol(x);norm(er,'inf')
%N=20;[u,x]=lap1mix17(f,a,b,al,dsol(b),N,mix=1,c);er=u-sol(x);norm(er,'inf')


 x=linspace(a,b,N+1)';
 h=(b-a)/N;
 A=sparse(diag(ones(N-1,1))-diag(ones(N-2,1),1));
 A=A+A'; %interior 3diag matrix with 2 on diag and -1 on sub- super diagonals
 F=(h*h)*f(x(2:N));
 %F(1)=F(1)+al;
 %F(N-1)=F(N-1)+be;
 F=[al;F;be];
 B=sparse(N+1,N+1);
 B(1,1)=B(N+1,N+1)=1; %mix =0
 switch(mix)
  case 0
     disp("Dirichlet");
  case -1 
     display("Left neumann; right:Dirichlet");
     B(1,2)=-1;F(1)*=-h;
  case 1
     display("Right neumann; left :Dirichlet");
     B(N+1,N)=-1;F(N+1)*=h;
 end    
 B(2,1)=B(N,N+1)=-1.0;
 B(2:N,2:N)=c*h*h*sparse(eye(N-1))+A;
 
 u=B\F; %all values
 %u=[al;u;be]; %we area adding bnd values

end