
function [u,x]=fdmbvpsolmix(f,a,b,N,c=0,al=0,be=0)
%a function solving BVP   -u''+cu=f in (a,b) u(a)=alpha u'(b)=beta
%by FDM on the uniform grid of N+! points in [a,b] i.e.
%we want to find [u(0),..,u(N)] s.t. u(k) \approx u(x(k)) for x(k)=a+k*h h=(b-a)/N
%INPUT
%f - a function handle to the rhs function f
%a,b- the interval endswitch
%N - no of gridpoints minus one
%c - coefficient of ODE (default =0)
%al be - boundary values at the ends (default =0)
%OUTPUT
%u - the solution vectorize
%x - vector of gridpoints
%EXAMPLE
%
%a=0;b=1;c=0;sol=@sin;dsol=@cos; f=@(x) (1+c)*sin(x);al=sol(a);be=dsol(b);N=20;[u,x]=fdmbvpsolmix(f,a,b,N,c,al,be);
%plot(x,u,x,sol(x));pause(1);plot(x,u-sol(x),";error;");
%sol=@(x) (b-x).*(b-x); f=@(x) -2*ones(size((x)); al=(b-a)*(b-a);be=0;N=20;[u,x]=fdmbvpsolmix(f,a,b,N,c,al,be);
%plot(x,u,x,sol(x));pause(1);plot(x,u-sol(x),";error;");
x=linspace(a,b,N+1)';
h=(b-a)/N;
C=A=sparse(N+1,N+1);
A(2:N,2:N)=sparse(diag(2*ones(N-1,1))) - sparse(diag(ones(N-2,1),-1))-sparse(diag(ones(N-2,1),1));
A(2,1)=A(N,N+1)=-1;
A*=(1/(h*h));
A(1,1)=1.0;
A(N+1,N+1)=1/h;
A(N+1,N)=-1/h;
C(2:N,2:N)=sparse(c*eye(N-1,N-1));

F=[al;f(x(2:N));be];
u=(A+C)\F;
end