
function [u,x]=FDmixlft(f,c,a,b,al,be,N)
#FD solver for -u''+cu=f u'(a)=al u(b)=be
#Input: f - function handle
#a,b- ends of interval
#N - no of mesh points -1 
#c - constant
#al, beta - Neumann and Dirichlet bnd conditions
#output:
#x- mesh points : (xk=a+k*h)
#u - FD solution on these mesh (interior mesh points - 3 point stencil) 
#f=@sin; u=@sin;a=0;b=1;
#[u,x]=FDsolmix(f,a,b,al=cos(a),be=sin(b),c=0,N=50)

if(c<0)
 printf("Warning: c=%g<0. The BVP -u''+cu=f may have no solution!\n",c);
endif
h=(b-a)/N;

x=a+h*(0:N)';
dd=-ones(N-1,3);dd(:,2)=-2*dd(:,2);
A=spdiags(dd,[-1,0,1],N-1,N-1);
A=A+h*h*c*speye(N-1);
#A- matrix for D bnd cond
L=spalloc(N,N,3*N-2);
L(2:N,2:N)=A;
L(1,1)=-h;
L(1,2)=h;
L(2,1)=-1;

F=zeros(N,1);
for k=2:N,
 F(k)=f(x(k));
endfor
F*=h*h;
F(1)+=al*h*h;
F(N)+=be;
u=L\F;
u=[u; be];
endfunction

